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ABSTRACT 



The Harris expansion method is applied to the elastic s-wave 
scattering of low energy electrons from hydrogen atoms and singly 
ionized helium atoms. The trial wave functions are Hylleraas 
functions of 22, 34, 50 and 70 parameters. It is found that rea- 
sonably accurate values of e-H phase shifts can be calculated but 
that e-He^ phase shifts are substantially less reliable. It is 
shown that the Harris method gives an accurate depiction of the 
location, but not the width, of the scattering resonances. Singlet 
and triplet s-wave phase shifts for e-H and e-He^ scattering are 
compared with the results of other calculations and H and He S 
state energy levels including the auto-ionizing levels are pre- 
sented and compared with other calculations and with experiment. 

It is tentatively concluded that the Harris method does not work 
on systems whose long range potential is of the Coulomb type. 
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INTRODUCTION 



I . 



The quantum theory of scattering has had a long and fruitful 
history, beginning with the very origins of quantum theory. Its 
extensive application to atomic collisions is indicated by the 
large number of texts, monographs and reviews of that subject pre- 
sently in print [l]. More recently the wide availability of high 
speed electronic computers has given additional impetus to the 
subject. New methods have been devised and older ones refurbished 
to take advantage of the computational power now at the disposal 
of virtually every investigator in the field. Parallel to the 
development of computer technology has been the application of 
powerful variational techniques to the scattering problem. 

Variational methods have held a distinguished place in mathe- 
matical physics since the days of the Bernoullis. Perhaps no 
technique has had wider application and thus it is no surprise that 
it has proven a powerful tool in scattering theory as well. Its 
first application to atomic scattering by Hulthen in 1944 [2], 
followed shortly by Kohn [3] ^nd Lippman and Schwinger [4] has led 
to extensive applications to all areas of scattering theory [ 5 l- 
One difficulty that arose with the Kohn and related methods how- 
ever, was the appearance of spurious singularities in the cal- 
culated phase shifts. While not fatal to the method they were 
annoying and some considerable effort has been expended in attempts 
to explain these singularities by Schwartz [ 6 ] and others [7,8,9]- 
By working to exploit rather than eliminate the singularities, 
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Harris [lO] developed an expansion technique that has yielded very 
reliable results to date. While this scheme involves some compu- 
tational simplifications over the Kohn method it too is not without 
disadvantages, which are dealt with in some detail below. 

Formal solutions abound to the scattering problem, ranging 
from the Born approximation and partial wave analysis to the 
Green’s function methods developed by Schwinger [ill, but accurate 
solutions to ’’real’' problems beyond the realm of potential scat- 
tering are relatively few, owing principally to the complexity of 
the problem when one permits the scattering center to have 
structure. Even the simplest of ’’non -tr i via 1” scattering problems, 
that of electrons or positrons incident elastically on ground 
state hydrogen atoms is already a ’’many body” problem and thus not 
subject to solution in closed form. 

It is the purpose of this thesis to investigate the application 
of the method developed by Harris to the problem of the elastic 
scattering of low energy electrons from hydrogen atoms and from 
singly ionized helium. Both have been attempted several times 
previously by various other methods, most notably in the case of 
hydrogen by Schwartz, whose definitive calculation was published 
in 1961 Cl2] and more recently by Callaway and his collaborators 
in 1969 and 1970 [l3jl4!l> and by Adelman and Reinhardt [l5l ^nd 
Truhlar and Smith [ 16 ] in 1972. Because of the more complex nature 
of the interaction in the case of singly ionized helium and a 
paucity of experimental data for comparison these calculations 
have not been as extensively pursued. The earliest attempt was 
by Bransden and E)algarno in 1953 Cl7] and the most fecent by Burke 
and Taylor in 1966 ClS] . 
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with the exception of Schwartz’ work roost of the important 
calculations of these problems have used the algebraic close cou- 
pling approximation C 19 ] wherein the trial wave function is made 
up at least in part of a combination of the first few bound state 
orbitals of the target atom. Such a scheme has earned a deserved 
reputation for accuracy and rapid convergence, however from a 
philosophical point of view it suffers because its application 
requires some a pr ior i knowledge of the detailed structure of the 
target particle. On the other hand Schwartz, in his calculation, 
made use of an ever increasing set of those functions first used 
by Hylleraas in his famous calculation of the ground state energy 
of helium r20] . This method makes no approximations other than 
those required by the use of a finite rather than an infinite basis 
space in which to work^. 

It was in the spirit of Schwartz’ work that the present project 
was undertaken. It was hoped that by utilizing the Hylleraas type 
functions in a calculation of electron-helium ion scattering that 
a more accurate and complete compilation of s-wave phase shifts 
would be obtained. As will be seen below such hopes were not 
realized, however it is hoped that the reporting of such results 



This is not strictly true. Electron capture by the scattering 
center is usually neglected in calculations of this type, and is 
neglected also in the calculation reported in this thesis. To 
incorporate capture requires including terms in the potential in- 
volving the interaction of the electrons with the electromagnetic 
field in order to account for the creation of the necessary photons 
and is properly a problem in quantum electrodynamics. The reported 
cross sections for capture in the system under consideration here 
are small compared to the elastic scattering cross sections and 
usually may be safely neglected. 



9 



I 



as were obtained will provide incentive for others to seek an 
explanation of the apparently anomalous results for helium ion 
scattering herein reported. 

In Section II a brief development of the partial wave method 
of scattering theory will be presented including the modifications 
required in the presence of the long range Coulomb force followed 
by a description of the Kohn method. Last will be a more detailed 
description of the Harris method and its application to elastic 
electron scattering by single-electron atoms and ions. In Section 
III the results of the numerical calculations are presented for 
both atomic hydrogen and singly ionized helium along with a dis- 
cussion of these results and comparison with those of other workers. 
An important by-product of the Fiarris method is the calculation of 
the bound state energy levels of the corresponding two-electron 
system. Results of these calculations for H and He are also 
presented in Section III. Most of the detailed calculations as 
well as discussion of the numerical methods employed are deferred 
to the appendices. 
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EXPANSION TECHNIQUES IN SCATTERING THEORY 



II . 

In this section will be presented a brief outline of the method 
of partial wave analysis [ 21 ] for the case of elastic atomic scat- 
tering. Part B will briefly present the Kohn method and in part C 
the Harris method will be developed in some detail along with its 
application to elastic He scattering of electrons. The system of 
units used throughout will be those in which m (electron mass), 
h (Planck*s constant) and e (electronic charge) all equal unity. 

In such a system the unit of energy is equal to twice the ground 

state energy of hydrogen (i.e., 27.2 eV) , and momenta are written 

2 

in terms of the wave number, k, such that E = k /2. It will be 
assumed that the nuclear mass, be it hydrogen or helium, is infi- 
nitely large compared to that of the electrons. 

A. PARTIAL WAVES AND PHASE SHIFTS 

— > 

Let a beam of particles of momentum k fall upon a spherically 

symmetric potential V(r). Taking the scattering center as the 

coordinate origin and the direction of k as the angular axis the 

wave function of the outgoing beam far from the coordinate center 

will be proportional to the sum of the wave function of the inci- 

— » — » 

dent beam, represented by the plane wave exp(ik-r) and a modified 
spherical wave, f (0) exp ( ikr )/r , representing the scattered portion 
of the incident beam. For proof that solutions to the Schroedinger 
equation of the form 

— ♦ — » 

vf / X ik • r ikr . , . . 

+ f(9)e A (1) 
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exist, the reader is directed to reference 21 (p. 220 ff ) . The 

differential scattering cross-section is defined to be 



0 (9)dA 



Number of particles per unit time scattered into 

the solid angle element dst at st 

Total Incident Flux 



( 2 ) 



where the particle flux is defined to be 



cp'' ^ - Yvy'') . (3) 

Putting (3) into (2) gives 

-4 ^ ^ 2 

cp • ds cp *r r 

' a(0)da = -S_ = _1_ Or. (4) 

o o 

- K ' 

where ds is the element of surface area at the detector and r is a 

unit vector in the direction of the scattered flux. cp and cp 

o s 

represent respectively the flux arising from the first and second 
terras of (1), respectively. Because of the scalar product in ( 4 ) 
only the radial part of the scattered wave function contributes 
to the scattered flux. Evaluating and cp^ from (3) ^nd applying 
them to (4) leads immediately to 



a(0)dA = I f (0) ! ^ dfl. . 



(5) 



It remains then to exhibit a means of calculating f(9) to 
complete the connection between theory and experiment. Such a 
means is the method of partial wave analysis. The most general 
solution to the Schroedinger equation in the case of axial sym- 
metry is 



Y(r) = E a, - Y (0) 

JL=0 ^ ^ ^ 



(6) 
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where is the spherical harmonic of order zero. The index i 

represents the angular momentum quantum number. Each radial function 
satisfies an equation of the form 




b / 2 b \ , , i(l+l) 

57 V rJ* ''I')* T2 

2r 



hf 1 

2 J r 



0 . 



( 7 ) 



If V(r) goes to zero fast enough at large r there v/ill be a region 
outside of which its effects can be neglected. Since is a 

function of r only, equation (7) can be rewritten 

u" + [k^-2V(r)- ^(f + l)/r^] u = 0 

and since the last two terms tend to zero for large r it is rea- 
sonable to assume that any solution u tends to the form A*sin(kr+c) 

i kr 

for large r. If this be so, u must have the form f(r)e where 
f(r) is approximately constant for large r, and (7) reduces to 

f" + 2ikf'- [2V(r)+ f (f+l)/r^] f = 0 

»l t 

and for large r, if f is nearly constant f « kf and this equation 
can be approximately integrated, giving 

2ik f = ^[2V(r)+ f (f +l)/r^]dr 

and the right hand side tends to a constant for large r only if 
V(r) goes to zero faster than l/r (ref. 22, p. 23). 

In the region where V(r) can be neglected equation ( 7 ) reduces 
to the spherical Bessel equation whose solutions are 

U. (r ) 
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where J^(kr) and ( kr ) are the spherical Bessel and Neumann 
functions, respectively C 23 ] whose asymptotic forms are 




sin(kr -i Tr/2)/kr 



^/2)/kr . 



2 

Since n^(kr) is not well behaved near the origin , if V(r) = 0 
everywhere then must be chosen to be zero in order for the 



solution u^ to be valid. Thus it can be inferred that the value 
of compared with will be a measure of the strength of the 
potential and hence of the scattering. To make this explicit (8) 
can be rewritten, using the asymptotic forms above 



so that equation (6) becomes for large r, a plane wave shifted in 
phase from that of an unperturbed (V=0) plane wave. can be 

taken equal to one and the normalization absorbed by the expansion 
coefficients in equation (6). 

The first terra on the right hand side of (1) can be expanded 
as follows: 




cos(kr- — ) j 



and letting B^/A^= - tan6 



^£( 07 i:^C^sin(kr- — + 6^) 



(9) 



00 



e 



ikr cos9 




(10) 



The Neumann function goes as 




b 
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and letting assume its asymptotic form, (6) and (1) can be 
equated and the result solved for f(0), giving, on choosing the a^ 
so as to assure an outgoing wave , 



f(9)=^E (2^+l)e 

^ £=0 



CO 



i6 



i 



sin6^P^(cos0) 



( 11 ) 



which is the usual form given for the scattering amplitude in terms 
of the phase shift. 

Deferring the calculation of 6^ to the next sections, the sub- 
ject of how the above formalism must be modified in the presence 
of a long range coulomb potential will now be considered. Recall 
that the classical solution for the motion of a particle in an 
inverse square central force field is a hyperbola. According to 
an argument due to Gordon (ref. 22, p. 54 ), if one considers the 
family of classical hyperbolic trajectories with one asymptote 
originating at the left limit of the z axis, the surface perpen- 
dicular to these hyperbolae goes as 



for large r. Thus it is as if the incident wave were initially 
distorted by the presence of the scattering center at infinity and 
its expected form would then be 



A similar argument with regard to the scattered wave leads to the 
form for the asymptotic wave function in the presence of a coulomb 



Z+ a k(r-z) = const 



e 




field 




ikz+ ia2m k(r-z) 
e ^ ^ 



+ f (Q)e : 



ikr- i^w': 2kr 



( 12 ) 



r 
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The coefficient a will be defined below. In solving equation (7) 
in the presence of the Coulomb potential one typically converts to 
parabolic cylindrical coordinates whereupon the Schroedinger 
equation becomes a version of the hyper geometr ic equation whose 
solutions are the hyper geometr ic functions (ref. 22, p. 57 ff ) . 

Using the asympototic forms of the solution functions as before, 
equation (9) becomes, in the case of the Coulomb potential 

fir 

— + a 07! 2 kr+ 77 ^+ 6 ^) (I 3 ) 

where 6^ is as defined above and 

= arg{r(f+l - ia ) } (l4) 

is the phase shift at infinity due to the Coulomb potential. 

Equation (11) must now be modified to account for this additional 
phase shift (ref. 22, p. 65 ff). 

The coefficient a in (12), (13) ^nd (l4) is proportional to 
the charge of the scattering center. However, since the problem 
at hand is one in which the pure Coulomb field exists only at long 
range, this charge is the net charge of the ion, in this case (Z-1) 
where Z is the atomic number of the scattering center, and in the 
system of units in use a can be written 

a = (Z-l)/k (15) 

and now 6^ can be interpreted as the phase shift due to the de- 
parture from pure coulomb scattering at close range. It is this 
quantity which is of physical interest. 

To summarize, it has been shown how the asymptotic form of the 
wave function (equations 1 or 12) can be calculated by solving the 
Schroedinger equation for a series of linearly independent functions 
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i-3 



(equations 9 or 13 ), in other words, that the problem can be reduced 
to solving a series of "partial wave" equations for the various 
values of the angular momentum quantum number A simple kine- 

matic argument indicates that as the energy of the incident beam 
is reduced the higher order terms in the expansion (6) become less 
important, enabling a fairly good representation of f(0) to be 
constructed using only the first few terras. 

B. THE KOHN VARIATIONAL METHOD 

Before taking up the Harris method it will be worth while to 
briefly examine the Kohn method since it is the most widely used 
variational method in low energy atomic scattering and its dif- 
ficulties were the motivation for Harris’ development. These dif- 
ficulties have boon examined elsewhere [ 6 - 9 ] and will only bo 
mentioned here. 

The Kohn method provides an estimate of the phase shift which 
is "second order accurate" in the error, i.o., the error term is at 
least second order small. In developing the Kohn method one assumes 
a trial function of the form 

CP = sin kr t cos kr + Y (i6) 



where t is an estimate of the tangent of the phase shift and Y is 
normally a linear combination of some basis set such that 



N 

Y = E a .X . -77^ 0 
i=l ^ ^ ^ 



( 17 ) 



so that cp is now a function of N+1 par ameter s - -the N values of a^, 
and t. Let be the true wave function whose asymptotic behavior is 






sin kr+ tan6 cos kr 



^(0)= cp(0) = 0 
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and define the error wave function 



e = cp - '3'. (iS) 

Now form the functional 

F = ^ cp(H-E)cp dr (19) 

which, using the fact that (H-E)'fcO, can be transformed into 

F = ^ cp(H-E)G dr ’ . 

Integrating twice by parts yields 



F = - - evcp) + 



dr 



^ e(H-E)cp 

and substituting (l8) for £ in the first two terms and for cp under 
the integral gives 

F = cp(H-E)c0 dr = |-(t-tan6)+ ^ e(H-E)C dr 



or 



k k 

- tan6 = - t- F + 






)e dr 



( 20 ) 



which identity is due to Kato [ 24 ] and exhibits the nature of the 
error terra. Now the Kohn prescription takes as its functional 



k k r 

I = — tan6 ~ 2 ^ ^ “ \ ^(H~K)cp dr 

and seeks its stationary value subject to the conditions 

bl bl 



( 21 ) 



bt ba 



= 0 , i=l ,N. 



The value of I thus found is the estimate of tan6 to second order 



C. THE HARRIS EXPANSION METHOD 
1 . Formal Solution 

It can be shown that the development of expressions 
for the phase shift lose no generality when the system, is restricted 
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to one of s-waves (^=0) scattered from a central potential. There- 
fore this restriction will be imposed on the present derivation of 
Harris’ method [lO] and the validity of its extension to atomic 
systems- will be assumed. 

With the above considerations in mind, assume a system 
describable by a potential V(r) and a Hamiltonian H. Construct a 
trial basis space using a set of N linearly independent functions 
In this restricted space the Hamiltonian can be represented 
by a matrix whose elements are 

H. j = <XilH|x^> (22) 



while the inner products of the basis vectors are 



L. . 



<xjx^> 



Then the eigenvalue equation 



(23) 



(H - Xl)C = 0 



where the elements of H and L are as given in (22) and (23) , can 
be solved for the N eigenvalues and corresponding eigenvectors 
C , (p =1 , . , jN ) . From the eigenvectors a new basis set is con- 
structed (which spans the same space as the Ix^;^) 



N 

|cp > = E c. lx 
^ i=l ^ 



(25) 



Now if the true wave function is represented by 

^ = S+tC + ^ 

where ^^(0)= 0 

^( ^ ) - s in kr + tan6 cos kr = S+tC 
r 
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(where now the identification of the terras S, t and C is clear) 
and if in addition | raay be reasonably well represented at sorae 
energy E by a linear corabination of the 

0)=1 ^ 



then the approxiraate wave function raay nearly satisfy the 
Schroedinger equation. If this is to be the case, then in the 
space spanned by the 'the vector 

(H-E) I S+ tC+ 



(26) 



raust be the null vector, which is, of course, merely the statement 
of Schr oedinger ’ s equation in the restricted space of the trial 
functions. This condition raay be satisfied if it is required that 
(26) have no component in the space spanned by the that is, 

that 

<cpp| (H-E)|s+ tc+ §> = 0, ( 27 ) 

Since the cp > are themselves linear combinations of the original 
0 

basis set 'the condition ( 27 ) is equivalent to requiring that 

(26) have no component in the space spanned by the Bearing 

in mind that the solution sought involves finding t, an estimate of 

tan6 , note that if <CP |H"e 1§'^ and <cp |h-e|s+ tC^ are separately 

r P 

equal to zero then (27) is satisfied and t can be immediately found. 

Note further that if E = X , the eigenvalue found in (24) corre- 

H 

spending to the vector that then 

r' 

<cp^| (H-E) U> = 0 ( 28 ) 
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for 



N 



<9 I (H-E) U> = E 
P co=l 



N N 

b E c. E <X. 
j=l 1=1 



(H-^p)lx/C.p 



and the third sura in this expression is raerely ( 24 ) in component 
notation, whence 
N 

E <xj (H-Xp) IXj> = 0 , j=l,...,N 



irrespective of the value 
follows immediately from 

tan6 = t = - 



of CO. Hence (28) 
( 27 ) that 

<CP^| (H-Ep) |s> 
<9p| (H-E^)|C>“ 



is satisfied and it 



(29) 



and the problem is formally solved. Finally as the size of the 
basis set is increased (i.e., approaches a complete set of 

quadrat ically integrablc functions) | should more closely 
approach |^>^and t should converge to the correct result. 

Kolker [ 25 ] has shown that the Harris method can be for- 
mulated as a variational principle and that under certain con- 
ditions it can be combined with the Kohn method to give a minimum 
principle. The advantage of the Harris method is principally com- 
putational, in that no integrals of the form <s| (H-E)|c^ need be 
evaluated. However the contribution of these integrals is needed 
to provide the second order correction to the estimate of tan6 in 
variational methods dependent upon the Kato identity (20), so it 
can be seen that the Harris method is necessarily less accurate 
than the Kohn method, as was found by Nesbet [8]. It remains true, 
however, that the correction term can be made arbitrarily small by 
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increasing the size of the trial basis set, for as this set 
approaches completeness, the error must approach zero. 

2 . A pplication to Atomic Scattering 

In applying the Harris formalism to the problem of atomic 
scattering the Hamiltonian is written as follows 



H 




2 

1 



1 

2 





(30) 



where the subscripts refer to the individual electrons, Z is the 
atomic number of the scattering atom and r^^ “ 

inter -electr on distance. The zero of total energy in this system 
occurs with all three particles at rest and separated from one 
another by an infinite distance. The trial basis set is taken as 
the Hyller aas -type functions [20] 





/ n £ 


1 n> 


0 

( ra - - 


Ks 










I 112 e . 


1 


2/ 



(31) 



where n, i and m are integers, chosen subject to the constraint 



n + £ + m^M, (M=l,2,...) (32) 

and O' is a variable parameter. 

This form of the wave function allows for electron exchange, 
i.e., the interchange of roles between bound and free during the 
scattering process, and enables the wave function to be properly 
symmetrized when the electron spins are ant i -par al lei (spin = 0, 
or ^singlet” state) or an ti -syrametr ized when the spins are parallel 
(spin = 1 , or "triplet’’ state) thus satisfying the requirements of 
the Pauli principle [ 26 , 27 ]. This is the only way the spin enters 
the problem and so once the symmetry of the wave function is prop- 
erly chosen no further concern need be taken with the spin 
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coordinates of the electrons. Inclusion of the r term in the 
potential and in the wave function takes account of the induced 
polarization of the atom arising from the repulsion felt by the 
bound electron as the incident electron approaches. Finally, the 
exponential term forces the functions to zero as r^ or r^ increases 
so that all boundary conditions required for these functions in the 
last section are met. As will be seen below the parameter (y. is 
varied by inspection so as to optimize the results obtained. 

The number of functions in the basis set is found by taking 
all possible combinations of n , f , m subject to the constraint (32) 
for a given M, and eliminating all terms which duplicate the 
function or require it to be identically zero irrespective of the 
choice of sign. The results for 1^ M ^ 8 are shown in Table I. 

M No. terms in basis set 

Singlet Triplet 



1 


3 


1 


2 


7 


3 


3 


13 


7 


4 


22 


13 


5 


34 


22 


6 


50 


34 


7 


70 


50 


8 


95 


70 



Table I. Basis set size for various M. 

The choice of coordinate system is of importance in sim- 
plifying the computational details. For the present computation 
the optimum choice seemed to be the one in which the principle 
coordinates are the radius vectors of the two electrons and the 
angle formed by the radius vectors. The six coordinates governing 
motion of the center of mass and orientation of the system in space 
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ar e 



cyclic, and the Hamiltonian becomes [28J 






^ 2 ^ _2 i i 

brj * " irl * ^2 ^'2 U? 4 ^®12 



s in0 



12 bo 



12 



_Z Z ^ 1 

111 
1 2 12 



(33) 



where r^^ = (r^ + i| - 2r^t2 

the volxune element for integrals in this system is 



'2 - ^'l"2 '=°=®12r by the law of cosines, and 



= STt^r^di^rldr^sinQj^dQj^ . 



(ih) 



The advantage of this system is that all of the integrals arising 
from forming the matrices H and L (equations 22-24) can be done 
exactly, while those arising from <cp [ ( H-E) | S + tC ^ can be done 
exactly in the case of hydrogen and reduce to single integrals 
which then can be done numerically in the case of helium ions. 

Because the electrons are identical and thus formally 
indistinguishable it is not necessary to integrate over the entire 
first quadrant of the sufficient to cover only 



the region this is true one must exami 



ne 



the integrals generated in forming and using (31) smd (33)- 

Carrying out the indicated operations gives rise to 52 terms for 

each H. . and four terms for each L. .. Each term consists of an 

ij 

integral of the form 

A(v,X ,u) = ^dr^dr^ e 



V-2 X-2 [i-2 



XXX 
1 2 12 



(35) 



or 



- w . r ~ v-2 X-2 U-2 

B(V,X,U-) = ^ ^ cosO 



1 2 12 



12 ' 



(36) 
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that every term involving A(V,\>0) and B(v,X,0) always has a 
coefficient of zero. Since the angular integration of B(v,A >2) 
vanishes, B(v,X>2) = 0, so using the recursion relations proved 
in Appendix B it is sufficient to explicity evaluate only the 
A(V,X,2) in order to find all the terms required to generate the 
H. , and L . . . 



Setting p = 2 in (35) ^nd performing the angular integration 



over the range 0 to n , gives, in terms of the integration region 
outlined in Figure 1 below and the volume element { 3 ^) 



which is the mathematical equivalent of integrating (37) in the 



particles possible because the electrons are identical, the in- 
tegrals (37) are sufficient to cover the entire quadrant. This 
argument is extended below to the integrals arising in 



il 






2 



(38) 



region r^ ^ r^ as can be seen immediately by making the variable 
change r^ in (38)- Thus as a result of the exchange of 



<9 I (H-E) I S+tC> . 



Equation (37) may be integrated directly, giving 




( 39 ) 
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Figure 1. Integration Region for A(V,X)2). 

2 

The factor l6 n appears in every integral; it can be dropped from 
the calculation since it will merely give rise to the same constant 
factor in both H and L and will not affect the solution of the 
problem. Because of the v^ay the recursion relations are structured 
the factor \/OL always occurs in the form 1/(T^ ^ ^ and thus can be 
factored out of the integration and incorporated into the coeffi- 
cient of each term in the expression for H. . and L. .. Since this 

ij ij 

is the only way the parameter OL enters these integrations, evalu- 
ating all the integrals need only be done once at the beginning of 
the calculation for a given basis and the values thus found can be 
used over again as OL is varied. 

Unfortunately a similar simplification is not available with 
the integrals from 1 (H-E)fs+tC^. However as will be seen only 

a relatively few integrals must be done numerically for each value 
of OL and k and so the time required to complete the numerical 
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integration (which dominated the time required to evaluate the 
helium ion phase shifts in all but the largest basis sets) was kept 
within reasonable bounds. Since the present calculation is con- 
cerned only with the s partial wave (f=0) further discussion will 
assume I has been set equal to zero in all equations of the form 
(9), (13) or (l4). Furthermore, since the proper form for hydrogen 
results if ( 13 ) is taken as the asymptotic form of the wave function 
and Z is sot equal to one in (15), the remainder of the development 
will bo concerned exclusively with establishing the Harris method 
in the presence of a long range Coulomb field. 

Bearing those comments in mind, note that as r goes to 
zero (13) does not behave as well as one would like. Because of 
the presence of the logarithmic term in the argument of the sine, 
as r approaches zero the interval between zeroes of the function 
approaches zero and u(r) does not approach a definite limit. To 
force u(r) to zero a factor is normally chosen which goes to zero 
fast enough as r goes to zero to override the effects of any 
singularity which may occur, and goes to unity as r gets large in 
order not to affect the asymptotic behavior of the function. For 

3 

the present calculation the factor was chosen to be ( 1 -exp ( -O^r /2 ) ) . 

This choice was made to insure that u(r) would go to zero fast 
enough that its oscillatory behavior near the origin would not affect 
the accuracy of the numerical integrations involving the asymptotic 
part of the wave functions. 

V'Jith these preliminaries aside, the form chosen for the 
asymptotic part of the wave function may now be presented: 



27 



s+tc = G 



^^1-e ^ Sn 2 ki^+ rj') 



+ taii6 cos^kr 2 + ” j~ ^ 



O' 



- ^2f ■ 2 " 

+ e (l-e 



(l-e ^ ^)^[sin(kr^+ 2m. 2kr^+ 77 ) 



tan6 cos^kr^+ ^ 2kr^+ 



(40) 



ivher e 



= ar 



^K-¥)}- 



(41) 



Evaluation of the complex Gamma functions is discussed in Appendix 

D. Nor mali;jat i on factors in (40) have been suppressed because they 

do not affect the calculation, cancelling when t is evaluated. 

Following the prescription outlined in section 1, the 

estimate of tan6 is evaluated from theratio of the integrals 

< cp l(H-E) I S ^ and < CD I ( H-E) I C ^ at the energy E found by solving 
Hr r 

(24). Since the functions (H-E)jx^^ have been previously evalu- 
ated in solving (22) and (23), considerable simplification is 
achieved by taking advantage of the fact that FI is a hetmitian 
operator and solving instead the equivalent integrals 

N 






(42) 



and 






(43) 



whence t is found from the negative ratio of (42) to (43) • Note 
again that 6 measures only the departure from pure Coulomb 
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scattering arising in the region close to the target where the 
three particles in the system must bo treated explicity. 

Carrying out the indicated algebra one finds that 
<s|(H-E)|x^ has 52 terms constructed of integrals of the form 



A^(V,X ,u)= 



-2 \-2 u-2 
r ^ r , ^ 

2 12 



, a, a 

e 



/ ^ \ 



X 



Z-1 



sin(kr 2 + - ^ 2kx^ + Tj) dr ^ 



(4M 



. r V-2 A-2U-2 “^^‘*' 2 ^ ^2 "2^1 ( 

,X ,H)= ^2 ^12 ^ ® 



O' \3 



As(V 



X 



Z-1 



sin(kr^+ ^ 2kr^+?7) 



B^(V,X ,^i)= 



-2 X-2 U -2 2 ^ 

^^2 ^12 ^ 



a, 0; 



^ ^ \ 3 

^1 ’2^2 ( "2^2) 

e \l-e / X 



sin(kr2+ 2kr^+ ?? )cos0^2‘^^i<^^2 



B2(V,X , 




V-2 

1 



X-2 

"2 ^ 



U-2 

12 



~{Z+ 

e 




O' 

2^1 




a 

2 




X 



sin(kr^+ ^ 2kr^^+ 7] )cos0^2^’^l*^^2 

Four similar integrals with cos(kr + ....) replacing the sine terms 
in (44) - (47) arise in evaluating < c|(H-E)|x^ and are designated 
> A^, and , respectively. Note that all these integrals 

are of the same form as those in Appendix B Section 1 and hence the 
same recursion relations may be applied. 
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AS with the integrals of equation ( 34 ) the evaluation of 
(44) - (47) starts with the case p. = 2 and as before, the angular 
integration yields B^(V,X,2) = B2(v,X,2) = 0 and introduces a 
factor 2 in the coefficient of Aj^(v,X,2) and ^2) . Also as 

with the A(v,X,2), terras of the form A^(v,X,2) and A2(v,X,2) 
always exist together and thus satisfy the requirements that allow 
the integration to proceed over ^ r^ only, as illustrated in 
Figure 1 . 

Starting first with A^(v,X,2), after angular integration 



Aj(vA,2)- 



Z-1 



sin(kr2+ ^ 



Now referring to Figure 2, note that the order of integration can 
be reversed by choosing the elements of area dr^dr^ as indicated 
while the sarae region is covered. With this change A^(V,X,2) 
becomes 

A^(v,X.2) = ^ . 



^0 



CO 

sin(kr2+ ^ 2 kt 79)^ ^ 



and now the integral over r^ can be done giving 

v+A-i “2^2 

• e \l-e 



OOrr2 , V 

32^ V : 



a 






J . C 

l(Z+%)^ *^0 



(V-l)I(Z+2 ) 



^" 2 " 2 



X sin(kr2 + ^ 
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( 48 ) 



Finally, making the variable change x = (Z + d)r, 






32n^Vl 



(2z+a) (Z+q:) 



v+X+l 






wher e 



ax 



£_\ 3 



V yfi 2(z+Q:)N . fkx ^ Z-1 

j(y)= \ e X (^1-e * -r 



^ 2kx \ , 

0n + n jdx 

Z+Q: W 



when Z =2 (He ), equation ( 49 ) must be evaluated numerically 
ever when Z = 1,J7 =0 and the integral reduces to a standard 
Csol . Numerical evaluation of ( 49 ) is discussed in Appendix 



( 49 ) 

How- 

form 

D. 




Figure 2. Interchanging the r^ and r^ Integration. 



Proceeding now to the evaluation of A2(V,X,2), note that 
the r^ integration may be done directly and that evaluating it at 
the upper limit (^^) gives the same function of r^ as did the 
evaluation of the r^ integration at the lower limit {^ 2 ^ 
evaluation of A^(V,X,2) but with v and X interchanged: 
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\ / a \3 a 

X 32tt^X ] (f 2 rC~ , (, ' 2 ’^lj ~ 2 ^ 1 V . z-io 

’^’2)= -^- ^ ^ r^sxn(kr^-.--P^ 2kr ^ 



A2(V 



+TJ) 



V /■ 2 \i 1 C~i " 2 ^ 1 ) V+X-i . Z-l„ .1 

'L \Tz^J -(nTyT \ "1 sin(kr^-.-^^2kr^+77)) 

i=0 ^ 



Making the variable change x = Oit^/2 in the first term and 
X = (Z+tt)t^ in the remaining terras, and recognizing that the 
summation is in fact A^(X,V,2) gives 



A,(V,X,2)= 



v+X+1 



( 2 z +a ) ^ ( 2 z +a ) 



j- A^(X,V,2) 



(50) 



wher G 



i.;(y) 



= fc"^ x>'(l-c-^)3sin( 



2kx 

"a 



. 2-^ ^ 



(51) 



Now since A^(v,X,2) and A2(X,V,2) always occur together it is con- 



venient to define 

Ag(X,V,2) = A^(V,X,2) +A2(X,V,2) (52) 

and similarly 

Ag(X,v,u) =Aj^(v,X,ii) +A 2 (X,v,m.) (53a) 

Bg(X,V,tl) = B^(V,X,^) + B2 (X,v,m.) (53b) 

A^{\,V,\i) =A^{v,X,u,) +A^(X,V,^0 (53c) 

B^(X,v,m.) = B2(v,X,M.) + B^(X,V,^X) (53d) 



where (53a) and (53b) arise in evaluating < s| ( H-E) |x ^ and (53c) and 

(53d) arise in evaluating <c|(H-E)|x^ and the number of terras in 

each is nov; reduced to 26. Unfortunately, explicit expressions for 

A, (v,X,2) and A.j(v,X,2) are required in order to evaluate A (v,X,l), 

■1-3 ® 
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B (V,X>1)> A (v,X,l) and B (v,Xj1)j respectively so the reduction 
s c c 

in labor is not so great as it might appear. Note that the factor 
2 

32tt /(2Z+Ck:) is common to every term, so it may be dropped from the 
calculation . 

3 . Summary of the Harris Method 

To summarize the Harris prescription as it applies to 
atomic scattering, the computational scheme is as follows: 

a. For the given basis set size evaluate the required 
"close-iny integrals from 



A(V,X ,2) 





(54) 



b. Evaluate the required A(v,XjM-) ^nd using the 

recursion relations of Appendix B. 

c. For a given O', evaluate H. . and L. . using the expressions 

xjxj 

. . . . V “^X 

in Appendix A. Note that each A and B must be divided by O' 

to give the correct result. 

d. Solve the eigenvalue problem (24) for the eigenvalues 
Ep and corresponding eigenvectors c^ . 

e. Find an appropriate incident electron momentum from 
the expression 



K = 2E^ 
P P 



+ Z 



(55) 



where -Z gives 
in rydbergs (1 
f . For 
of a, evaluate 
required values 



the ground state energy of the scattering center 
ryd = 13.6 eV) . 

this value of k and the previously assigned value 

(numerically if necessary) l^, and for the 

of the index y using (49) ^nd (51) for I and I and 

1 s 
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1 



I 



'• M Z 



eo Q:x 

-X y{^ ” 2(z+0i) Y ( kx Z-1 ,, 2kx \ 

x'll-e ' ’ ) cos [ — — + — — 071 + n) 

, \ / Vz+a k Z+Ck; / 



,(y)= ^ x^( 



dx ( 56 ) 



and 



^ r “X y,_ -x*3 /2kx Z-1 n 4kx 

x'^(l-e ) cos ^ “ + T]jdx 



(57) 



g. Using these values, evaluate the required terms of , 

^1’ ^3 ^ = 2 using (48) and (50). 

h. Using the recursion relations in Appendix B evaluate 

the remaining required terms of A , B , A and B . 

^ s s c c 

i. Evaluate the < sI(H-E)!x-^ and < cl(H-E)|x-> using 

N 

the expression in Appendix A and form E | ( H-E) | X ^ ^ ^nd 

c <s| (H-E)lx. > . 

j=l IP 1 

j. Finally the negative quotient of the two sums above is 



tan6(k ), the desired result, and repetition of the above step for 
r 

various OL normally permits the entire desired range of k^ to be 
cover ed. 
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III. RESULTS AND DISCUSSION 



The calculations reported in this section were performed in 
double-precision ( l 6 figures) arithmetic on the IBM 36O/67 computer 
of the Naval Postgraduate School Computer Facility. The programming 
language v;as Fortran IV and the H-level compiler was used. The sub- 
routines used to factor the matrix and invert U and to solve the 
eigenvalue problem, and the 32 -point Gauss -Legendr e quadrature 
subroutine, as well as several auxiliary subroutines were furnished 
by the source library of the Computer Facility. 

The principal program was, of course, written to evaluate the 
phase shifts. This program and its required subroutines are re- 
produced following the aj)pendices . In addition it was found con- 
venient to have available modifications of the main program which 
solved only the eigenvalue problem, giving in one case the values 
of the incident electron momenta in the elastic scattering range 
and in the other the calculated energy levels of the bound state 
system. Since the scheme adopted for checking the eigenvalues and 
eigenvectors was rather time consuming and involved a rather large 
quantity of printed output, it was done on only selected matrices 
and then only in the supplementary programs. The modification of 
the phase shift program to perform these supplementary functions is 
obvious and the programs so modified are not reproduced here- 
Several features were built into the programs to permit the maximum 
flexibility with the minimum of internal changes and as it finally 
evolved only the storage capacity within the program had to be 
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changed, depending on the size of the basis set to be used. All 

other options were controlled by data inputs from external devices. 

The results for e-H scattering will be presented first, followed 

+ 

in part B by the e-He results. Comparison with the results of 
other calculations will be made where appropriate. The results 
will be presented in the following order: (i) phase shifts for 

singlet and triplet scattering; (ii) location of resonance levels; 

and (iii) bound state energy levels from the eigenvalue problem 
(see Appendix C) . 

A. ELECTRON -HYDROGEN SCATTERING 
1 . Phase Shifts 

Singlet and triplet phase shifts were calculated for basis 
sets from N = 22 elements (M = 4, singlet; M = 5, triplet) to 
N = 70 elements (M = 7, singlet; M = 8, triplet). Although the 
program had built into it the capability of calculating phase 
shifts for the 95“eleraent basis set (M = 8, singlet) the eigen- 
value problem began to show signs of instability at this size and 
the computer time required to find eigenvalues and eigenvectors 
became unacceptably long (8-10 minutes per matrix), permitting no 
more than testing the program at this level. In addition the 
storage area required for such programs approached the limit avail- 
able on the computer . 

Preliminary to calculating the phase shifts, equation (24) 
u’as solved for a wide range of values of the non-linear parameter OL 
and the resulting values of incident electron wave number were 
plotted as a function of O'. As expected, the eigenvalues fell into 
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clearly discernable trajectories whose evolution as N was increased 
was easily traceable, at least in those regions where the eigen- 
values were sparse. Figure 3 shows such a plot for the singlet 
states of hydrogen with N = 70 . 

This plot clearly illustrates a major disadvantage of the 
Harris method, which has been alluded to by others [l3>25l in the 
context of the inconvenience of being unable to pick the scattering 
energy without the expenditure of considerable labor. The dis- 
advantage seems to be of a more fundamental nature, however. While 
coverage of the entire energy range is certainly possible, at least 
in most cases (it was not possible with the basis sets used to reach 
zero energy in the e-H triplet case), sufficient points at any given 
value of k are not always available to permit the sort of investi- 
gation of convergence rate that characterized the analyses of 
Schwartz [l2j and Armstead [ 31 ]- Furthermore in those regions 
where few eigenvalues exist there is no guarantee that they will 
exist for an optimum value of O', Schwartz has indicated that the 
useable range of CC was between about 0.8 and 4.0 in his calculation, 
and this generally proved true in the present case also. Reference 
to Figure 3 shows that even in the largest basis set used the value 
of a. for the only trajectory passing k = 0.1 is around 5.0, just 
outside this range. In this region, increasing N is not likely to 
produce additional eigenvalues at which to solve the phase shift 
problem, for as is shov;n in Appendix C, the eigenvalues result from 
a Ray leigh-Ritz variational calculation of the energy levels on the 
two-electron bound system (see part 3 of this section) and hence 
the eigenvalue trajectories shown in Figure 3 represent approximations 
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Figure 3- Eigenvalue Trajectories for Singlet e-H System N = 70 . 
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to the various energy levels of H for a given (X. Although the 
scale is not the energy scale one usually associates with such 
plots the relationship between the two is 

= 2E + 1 (58) 

in the case of hydrogen, where E is given in atomic units (27.2 eV) . 
In such calculations, as N increases the lowest lying trajectories 
become quite good approximations to the energies of the lowest lying 
eigenstates of the bound system [ 32 ] after which further increase 
in N will' give rise to no additional trajectories in that region. 

In fact the situation is likely to get worse from the point of view^ 
of the scattering calculation, since as N is increased the minima 
in the trajectories become broader, approaching horizontal straight 
lines as completeness is approached. Consider an incident electron 
wave number near but above that corresponding to the highest well 
defined energy level of the bound system for a given N. As N is 
further increased the values of (X for which such a wave number is 
an eigenstate of the system will be pushed closer to, and perhaps 
beyond, the limits of the lange in CX for which reliable results can 
be expected. Figure 4 illustrates this behavior for a typical 
eigenvalue trajectory. As N is increased the curve becomes broader 
and its minimum approaches a limit. For N ^ 70 the minimum probably 
will not decrease appreciably but the curve will become broader. 

One might conclude that in an energy region which is rich 
in eigenvalues the situation ought to improve. That such is not 
always the case can be seen by examining Figures 5 and 6, showing 
singlet phase shifts for k = 0.4 and 0.8 plotted as a function of (X . 
It must be noted that the lines connecting the points arising from 
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Figure 4* Typical Evolution of One Eigenvalue Trajectory as N 
is Increased. 
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the calculations for a' given basis set are just that and no more. 

They have no physical or mathematical moaning, for the points 
plotted represent a complete collection of the possible values 
to be found at a given energy (possibly neglecting one or two 
at small Of) . 

In Figure 5, while the ’^curves’^ of all basis sets have 
clear cut minima, there is no clearly defined trend as the basis 
set size increases, as Schwartz and Armstead observed in their Kohn 
calculations. This may bo duo to numerical errors which accumulate 
more rapidly as the basis set is increased in size, however as will 
be seen below the ambiguity is considerably reduced in the triplet 
case ivhich loads one to conclude that this error is probably small. 
Note however that the range of Of over which the phase shift is 
nearly constant is greatest for the case N = 70 and hence the 
minimum of that curve has been selected as the most probable value 
in this calculation. The rather large error indications in Table 
II reflect this situation. 

Referring now to Figure 6, note that while the minimum of 
the N = 70 curve is at 6 = 0.761, there is a plateau forming 
around 6 = 0.886. For this value of k increasing N may well bring 
improved results since the minima of more trajectories may be ex- 
pected to move into the region, however the ambiguity will remain 
if the plateau becomes broader and stable in value. A similar 
structure is evident in the calculations for k = 0.7 smd 0.6 as well. 
It is interesting to note that these plateaus correspond to the 
values for the phase shift given by Schwartz and the difference 
between these plateaus and the minima of the curves increases with 
increasing k until at k = 0.866 the difference is in excess of 10%. 



6 ( radians) 




a 



Figure 5 - Singlet Phase Shift for e-H as a Function of (X 
k = 0.4. 
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Although they do not discuss it, close examination of the results 
of the Harris calculation of e-H singlet phase shifts by Oberoi 
and Callaway C 13 ] using the close-coupling approximation seems to 
indicate that their results may be as much as 4-5% less than 
Schwartz’s at k = 0.8, although their graphical presentation of 
results makes such a comparison difficult. 

Appeal to experiment to resolve the ambiguity is fruitless 
at this time for the difference in total scattering cross-section 
at k = 0.8 for the two values of 6, 0.886 and O .761 is only about 
3 %, far less than the accuracy of any experiments conducted to 
date [ 33 ]- The experiments report results about 15% less than 
those of Schwartz and Armstead at the upper energy limit and it is 
worth noting that the present result is also lower than that of 
Schwartz’s. However because of the large experimental errors re- 
ported 20%) the significance of this fact, if any, cannot be 
estimated . 

Based on the results available it seems more wishful than 
logical to ascribe more significance to the plateaus in the N = 70 
curves than to the minima and in the absence of the results of 
Schwartz the lower value would most likely have been chosen as the 
most probable value. Therefore the minima will be chosen for 
cons is tency . 

The extrema of the curves for phase shift vs. Oi generally 
are minima except at the lowest energies calculated and they seem 
to trend dov;mvard, so in most cases the value of the phase shift 
for N = 70 can probably be taken as an upper limit except for 
k ^ 0.2 where it seems to be a lower limit. Where an unambiguous 
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Figur e 6 . 



Singlet Phase Shift for e-H as a Function of Oi 

k = 0.8. 
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convergence could be established it seemed to be about as the 
sequence 1/N. 



The scattering length is defined to be 

tan6 

o 

a = lim — r 

O K 

k-*0 



( 59 ) 



In the present calculation it seemed natural to evaluate the 
scattering length for e-H scattering by simply evaluating (tan5)/k 
for a succession of small values of k and extrapolating the results 
to k =0. This was done for the case N = 70 and the smallest value 
of k used was .0004l2. The extrapolation to k = 0 gave 

a = - 6 . 02181 ( 1 ) 

where the figure in parentheses is the uncertainty in the last 
digit. This value differs by about 1% from the value = -5.965 
reported by Schwartz [l2j. Note that the uncertainty reported 
above is the uncertainty arising from the extrapolation to zero 
and is not necessarily an indication of absolute accuracy. The 
extrapolation used in evaluating a^ is shown in Figure 7. 

It should be noted that the difficulties with the 
Harris method discussed above should not be j^eculiar to this choice 
of trial wave function, for since the Hamiltonian (30) is the exact 
Hamiltonian for any two-electron system the eigenvalue calculation 
should lead to the same effects described above for any choice of 
trial wave function as long as it meets the usual boundary condi- 
tions and is based on a quadra tica lly integrable basis. Perhaps 
by using some approximate Hamiltonian the eigenvalue situation 
would be improved but the loss of accuracy from this approximation 
would probably offset any gains from improved eigenvalue behavior. 
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Evaluation of e-Fi Singlet Scattering Length. 
N = 70 
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The results for the triplet phase shifts are much as has 

been described above for the singlet case except that the curves 

seem to trend upward toward maxima for values of k larger than was 

found in the singlet case. Since in the triplet case the electron 

spins are parallel the effective repulsion between the electrons 
due to the Pauli principle tends to improve the convergence. Thus, 
while the qualitative picture is similar to that of the singlet 
scattering, the ambiguities are much less severe. Figure 8 shows 
the eigenvalue trajectories for N = 70. Note that there is no tra- 
jectory that reaches k = 0, the lowest one just barely passing 
k = 0,1, Hence it was not possible to evaluate the triplet scat- 
tering length as was done for the singlet case. Figures 9 ^md 10 
show plots of phase shift vs. OL for k = 0.4 2 md 0.8. Finally, the 
curve of 6 vs . k for both singlet and triplet scattering is shown 
in Figure 11. The resonances shown on the plot are discussed in 
the next section. In each case one resonance which has been re- 
ported elsewhere as being just below the excitation threshold was 
observed in this calculation to be just above and is so indicated 
on the plot by the vertical dashed lines. It is believed that in- 
creasing the basis set size would bring these values below the 
threshold. The singlet and triplet phase shifts are tabulated and 
compared with those of Schwartz [ 12 ] in Table II. 

2 . Scatter i nq Resonances 

The present calculation offers a feature not found in 
Schwartz ^s calculat ion--the localization of real scattering re- 
sonances. There is probably no theoretical reason why the Kohn 
method should not show the resonances, however since they are very 



47 




Figure 8. Eigenvalue Trajectories 
System. N=70 



for Triplet e-H 
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Figure 9. Triplet Phase Shift for e-H as a Function of O'. 

k = 0.4 
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Figure 10 



Triplet Phase Shift for e-H as a Function of a. 

k = 0.8 
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Figure 11. Calculated Phase Shifts for e-H Scattering. 
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narrow [i 8 , 34 j 353 it is quite likely that they simply were not 



resolved. It 


might also 


be difficult to 


distinguish them 


from the 


spurious ones 


characteristic of such calculations. 




Incident 










Electr on 




Scattering Phase 


Shift in Radians 


Wave Number 


Singlet 


Tr iplet 






Pr esent 


Schwar tz[ 12] 


Pr esent 


Schwar tz[ 12] 




Work 


Wor k 


0,1 


2.566(5)^ 


2.553 


2.9386(5) 


2.9388 


0.2 


2.067(4) 


2.0673 


2.7187(7) 


2.7171 


0.3 


1 .696 (6) 


1 .6964 


2.5013(7) 


2.4496 


0,4 


1.414(6) 


1.4146 


2.2980(5) 


2.2938 


0.5 


1.194(7) 


1.202 


2.1093(3) 


2 , 1046 


0.6 


1 .020(4) 


1.041 


1.9321(5) 


1.9329 


0.7 


.878(3) 


.930 


1.7772(4) 


1.7797 


0,8 


.761(4) 


.886 


1.6389(5) 


1.643 


0.9 


.693(7) 


— 


1.5622(6) 


1.563 



^The figures in parentheses represent the uncertainty in the last 
digit reported. 

Table II. Singlet and Triplet Phase Shifts for e-H 
Elastic Scatter ing . 

The existence of scattering resonances has long been 
associated with the auto-i onij'ing states of the compound system [36] 
and recent experiments and calculations have confirmed this 
l' 35, 37,38]. Aut o-ionizing states (in the case of H they are more 
properly called auto-detaching) are those in which the total energy 
of the system exceeds the ground state energy by more than the 
ionization energy of the two-electron system but in which both 
electrons are in excited states. Such states can decay non- 
radiatively when one electron's excess energy is transferred to 
the other which is then ejected from the atom while the first 
reverts to the one-electron ground state. Such states are normally 
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extremely short lived, typical lifetimes often being several orders 
of magnitude less than the normally allowed radiative transitions C39] . 
Clearly an electron captured momentarily into such a state will be 
considerably delayed in its passage by the atom, suffering a cor- 
respondingly greater phase shift, and thus show a resonant behavior 
in the cross section. Since these states lie above the one electron 
ionization potential and are extremely short-lived they are often 
referred to as quasi-bound [4o]. 

Figures 3 ^md 8 show the evidence of such states in the 
singlet and triplet structure of H . They are the broad shallowly 
curved trajectories lying just below the excitation threshold in 
each case. 

Since these states do lie above the ionization potential 
of H and thus among a continuum of states, one expects from the 
Hylleraas -Undheim theorem [ 32 ] that such states would not be well 
approximated by an eigenvalue trajectory until after all lower 
lying states had been. The reasons for this rot being so are not 
well understood, however Holr^ien and Midtal [4o] have found that if 
the trial wave functions include an appropriate mix of the one 
electron state functions (i.e., Isns , n > 1, and also Isis, 2s2s , 

2p2p, etc.) that the eigenvalues corresponding to the autoionizing 
states begin stabilizing long before those of many lower lying 
states. They conjecture that this probably occurs because the 
equal treatment of the two electrons leads to an approximate 
or thogonalization to the lower states. The limits approached by 
these trajectories may not be upper bounds however. 
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That such conditions occur in the present problem is most 
likely a result of the use of Hyller aas -type trial functions, which 
can be combined to approximate the one-electron orbitals needed to 
satisfy the condition set by Hol/Sien and Midtal above. 

In the present calculation two resonances were resolved 
within the e-H elastic scattering range. One in the singlet states 
was observed at k = 0,83713 which corresponds to a bound state 
energy of E = -.l 496 lO atomic units (1 a.u. = 27.2 eV) . A second 
in the triplet states was observed at k = 0 . 864 l 8 which corresponds 
to E = -.126596 a.u. In addition there were two such states identi 
fied just above the excitation threshold which probably correspond 
to levels predicted just below the threshold by others [38]. In 
all liklihood a larger basis than those used here would give a 
better approach to the results reported by others. Since the means 
of locating these resonances was from a bound state calculation no 
estimate of their width can be given. A characteristic of these 
calculations is that the phase shifts associated with a particular 
eigenvalue trajectory tend to be rather independent. Thus evaluati 
a phase shift at a value of k corresponding to a resonance energy 
but at a value of O' corresponding to a normal (i.e., non -r es onant ) 
eigenvalue trajectory appears to show little or no evidence of the 
resonance. Although the non-resonant and resonant trajectories 
sometimes do cross, no conscious attempt to evaluate the phase 
shift at such a crossing point was made since in the eigenvalue 
method used (Jacobi variable threshold C 4 l] the eigenvectors found 
for degenerate or nearly degenerate eigenvalues may be grossly in- 
accurate (Ref. 4 l, p. 281) and thus the phase shift calculated at 
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this point would be difficult to interpret. The independence of 
the phase shifts corresponding to different trajectories is probably 
related to the orthogonality properties inherent in the eigenvalue 
problem . and discussed above but beyond this it is not well understood. 
Table III summarizes the results for the e-H resonances and compares 
them with other calculations and with such experimental data as is 
available . 

3 * H Energy Levels 

As shown in Appendix C and as demonstrated in the last 
section, the first part of the Harris pr escr ip>t ion involves what 
amounts to a variational calculation of the bound state energies. 
Reports of such calculations abound in the literature, commencing 
with Hylleraas [20] and continuing to the present, and have 
attained a high degree of accuracy. Hence the results obtained 
here are of no more than passing interest except that they arise as 
a natural by-product of the Harris method, and this fact seems not 
to have been mentioned previously in the literature concerning the 
Harris method. Considering that no special pains were taken to 
assure accuracy (other than the normal and reasonable ones dis- 
cussed in Appendix D) the energy results are surprisingly accurate, 
and since they are obtained with no additional effort from the 
phase shift calculation it seems worthwhile to discuss them briefly 
her e . 

Peker is C46] has reported a calculation of the ground state 
energies of He and H in which he used matrices of up to order IO 78 
and was at great pains to achieve the maximum possible accuracy. 

He reports the ground state of H from this calculation to be 
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-.52775097 a.u. In the present work a value of -.5277477 a.u. was 
obtained with an order 70 matrix. 

As an example of the bound state results obtainable, 

Figures 12 and I 3 show the negative energy structure for the 
singlet {N = 70) and triplet (N = 50) calculations. Not all the 
detail is included in the region near zero energy because the large 
number of eigenvalues in that region makes it very difficult to 
sort out the eigenvalue trajectories. Note that the lower eigen- 
values clearly exhibit the broad minima alluded to in section 1. 

f 

This characteristic of the trajectories is even more pronounced in 
the case of helium, which will now be discussed. 

B. ELECTRON -HELIUM ION SCATTERING 
1 . P hase Shifts 

In principle the only difference between scattering elec- 
trons from helium ions and from hydrogen is the presence of the 
long range Coulomb force in the former. This effect is accounted 
for in the form chosen for the asymptotic part of the trial wave 
function and once the more difficult problem of numerically evalu- 
ating the resultant integrals is solved, one expects the procedure 
to be straightforward and to yield results whose accuracy is de- 
graded only by the relative inaccuracy of numerical integration vs. 
exact integration. Appendix D contains extensive discussion of the 
method of evaluating the numerical integrals and assuring their 
accuracy, and it appears that a high degree of confidence can be 
placed in them. Nevertheless the results of these phase shift 
calculations are not so accurate as the hydrogen ones, if one 
accepts the results of Burke and Taylor Cl8] as the most accurate 
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Figure 12. Singlet Energy Levels for H~. N = 70 
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Figure 13- Triplet Energy Levels for H . N = 5^ 



59 



to date. The close agreement of their results with those of 
others by different methods [35,44] indicates that this confidence 
is probably well placed. Thus it would appear that the Coulomb 
force has an effect on the accuracy of the calculation beyond what 
is initially expected. 

As Figures l4 and 17 show, the eigenvalue trajectory situ- 
ation is somewhat improved in the case of helium, owing to the 
existence of a complete set of two electron bound states in the 
region below the first ionisation potential of He. Nevertheless 
the singlet phase shifts at low energy appear to be as much as an 
order of magnitude larger than the results reported by Burke and 
Taylor . Because the calculation gives only tan6 the value of 5 
is uncertain by ^ n^, and the choice of 6 to be .4 at k = 0.2 is in 
closer agreement with Burke and Taylor. However the slope of the 
6 vs . k curve at k =0.2 was measured and its value seems to support 
the higher value chosen. At higher energy the convergence appears 
to be to a result about 5-10% higher than that reported by Burke 
and Taylor, although the existence of the scattering resonances 
at higher energies undoubtedly affects the results and makes com- 
parison difficult in this region. The existence of the resonances 
makes more likely the possibility of attempting to evaluate the 
phase shift at a nearly degenerate eigenvalue, with its concomitant 
inaccurate eigenvector which may further confuse the picture. 

Typical plots of 6 vs . O; for singlet e-He^ scattering are sho^vn in 
Figures 15 and l6 . 

A similar but (not unexpectedly) less severe situation 
exists with the triplet scattering. Figures l8 and 19 show typical 
triplet plots of 6 vs. O: . 
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Figure l4 . 
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Eigenvalue Trajectories for Singlet e-He 
System. N=70 
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Figure 15 . Singlet Phase Shift for e-Hc^ as a Function of O', 
k = 0.8 
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k = 1.4 
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Figure 17. 



Eigenvalue Trajectories for Triplet e-He^ 
System. N=50 
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6 ( radians ) 




Figure l8. Triplet Phase Shift for e-He^ as a Function of Ci . 
k = 0.8 



65 



6 ( radians ) 




a 

Figure 19. Triplet Phase Shift for e-He as a Function of O'. 

k = 1 .4 
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Results of the present calculation are presented and com- 
pared with those of Burke and Taylor in Table IV and plotted in 
Figure 20. 

2 . .S catter inq Resonances 

It was possible to identify six singlet resonances and four 
triplet resonances. As shown in Table V the agreement with the re- 
sults of Burke and McVicar [44l and others [37,40,42,43)4?] is 
quite good. 

Burke and McVicar give the energies of seven and five re- 
sonances for singlet and triplet states, r espectively . As in the 
hydrogen case an additional resonance level could be identified just 
above the excitation threshold by examining the energy level plots 
for each system. Again as in the hydrogen, these levels would 
probably appear just below the excitation threshold if the calcu- 
lation were done with a sufficiently large basis set. 

3 . Helium Energy Levels 

The helium energy spectrum has been even more extensively 

calculated than has hydrogen [20,46,48,491. The most accurate such 

calculation is that of Pekeris [46] who has reported a value of 

- 2.903724375 a.u. for the ground state energy of He. After applying 

mass polarisation and relativistic corrections and correcting for 

< - 1 

the Lamb shift, he finds a value of I. 983 IO 687 x 10 cm for the 
ground state of He, which compares with the experimental deter- 
mination by Herzberg [ 50 ] of I. 983 IO 8 x 10^ ^ .05 cm Similarly, 

for the 2 S state of He, Pekeris obtained -2.17522937822 a.u. or 
3.84566 X lO^cm ^ while Herzberg measured 3-845473 x 10^ + .05 cm 
In the present calculation -2.90372410 a.u. and -2.175191 a.u. were 
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Scattering Phase Shift in R adians 



Incident 
ElGctr on 
V^ve Number 

0.2 

O.z* 

O.A472 

0.6 

0.6325 

0.7746 

0.8 

0 . 8944 

1.0 

1.1832 

1.2 

1.3416 

1 .4 

1.4832 

1.5811 

1.6 

1.6553 
1 .6712 
1.6971 
1.7073 
1.732 



Singlet 
Pr esent 
V^or k 

3.5(9)^ 
.98(4) 

.63(2) 

. 512 ( 8 ) 
.446(7) 
.417(7) 
.407(4) 

.348(4) 

.326(3) 



Burke and 
Taylor [l8] 



.4217 

.4092 

.3989 

.3912 

.3850 

.3780 

.3779 

.3910 

.2977 

.3988 



.2572 



Tr iplet 
Pr esent 
Wor k 

4.1(9) 

1.469(8) 

1.108(6) 

.965(6) 

.889(4) 

.831(5) 

.778(4) 

.730(4) 

.705(4) 



Burke and 
Taylor [18J 



.9088 

.8873 

.8649 

.8462 

.8253 

.7895 

.7593 

.7341 



.7006 

.6931 



^The figures in parentheses represent the uncertainty in the last 
digit reported. 



Table IV. 



Singlet and Triplet Phase Shifts for e-He 
Elastic Scattering 
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Figure 20. Calculated Phase Shifts for e-He^ Scattering. 
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“ionizing Levels Corresponding to e-He Scattering Resonances. 




m. 











obtained 
21 and 22 
triplet S 



13 

for the 1 S and 2 S levels in He, respectively. Figures 
show the calculated energy levels for the singlet and 
states of He. 
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Figure 21. Singlet Energy Levels of He. N=70 
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Figure 22. Triplet Energy Levels of He. N=50 
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IV. CONCLUSION 



The principle advantage in calculating phase shifts by the 
Harris method is that it is not required to evaluate matrix 
elements involving only the asymptotic part of the wave function. 
This results in a considerable computational simplification but 
the price paid is that the error is now of first order and hence 
the method is inherently less accurate than methods dependent on 
the Kato identity. 

An additional disadvantage is that in certain regions there may 
be a paucity of energy eigenvalues at which to carry out the calcu- 
lation. Since these eigenvalues are approximations to the compound 
system energy levels, in those regions where the approximations are 
good, increasing the number of parameters in the trial wave 
function may simply make the situation worse. Even in regions rich 
in eigenvalues the fact that they are related to the bound states 
may have unpredictable effects on the phase shift calculations. 

The Harris method does have a feature that could prove useful 
in complex systems. Resonances in the scattering cross-section 
can be uniquely identified as arising from certain bound states of 
the compound system, such as the au t o- i on iz ing states in the pre- 
sent instance. Since the scattering at a particular energy can be 
identified v/ith the compound state energy level, when that level is 
well defined by the trial wave function the correlation with the 
behavior of the phase shifts in that region will be clear. 
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+ 

\ChGthor the anomalous results for the e-He system found here 
are an indication of a defect in the Harris method or due to a 
subtle error in computation or analysis on the part of the author 
is not absolutely clear. Although the latter cannot be ruled out, 
the success of the method in calculating the e-H scattering using 
the same computer program and even, as pointed out in Appendix D, 
the same integration scheme (on a test basis), seems strong evi- 
dence in favor of the former conclusions. A third possibility, of 
course, is that the results obtained here are correct and the re- 

I 

suits of Burke and Taylor and others are incorrect. In the 
absence of experimental evidence to the contrary, the consistency 
of the other calculations would seem to make this conclusion also 
unlikely . 

Thus it is the tentative conclusion of this thesis, that the 
Harris method should be applied to electron-ion scattering with 
some degree of caution. 
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APPENDIX A 



Formulae for the Various Matrix Elements 



1 . H . . and S . . 

U 

The matrix elements < x.|h1x*^^^ made up of 52 terms, each of 



which includes an integral of the form 



, r - - ■ V-2 X-2 H-2 

A(v,X .Li)= 



(A-1) 



or 



V . r v-2 X-2 [1-2 

B(V,X ,ii)= r 



12 cos9 



12 ■ 



(A-2) 



The trial wave function is 



U 



X = G 






/ n f f n\ 
V1^2 - ^1^2A 



m 

12 



(A-3) 



and can be characterized by the three exjjlonents n,f ,m. Let )( ^ be 
character ized by n . , f . , m . and X • t>y n . , I . . m . . With this 

1 1 1 j j j 

notation H. . can be written 
ij 

H. . = - (0'^/4) •Afn.-*- n.+ 2,f.+ f.+ 2,m.+ m.+2) 

J J ' ^ X J X J ' X J ^ 

+ [ ( n . + 1 ) -Z ] *A(n.+ n.+ l,i'.+ 2,m.+ m.+2) 

J 

+ [^(f .+ 1) -Z ] 'Afn.^- n .+ 2,f.+ f .+ 1 m.+ m.+2) 

+ ^m.*A(n.+ n .+ 2,f.+ f .+ 3,ni.+ m. ) 

3 1 J ’ 1 J ’ 1 J ' 

+ ^ra.*A(n.+ n .+ 2,m.+ m. ) 

J 1 J 1 J 1 J 

- m.(m.+ f.+ n.+ l)*A(n.+ n.+ 2,f.+ f.+ 2,m.+ m. ) 

J J J J ' ^ X J ’l J ’i J ^ 



%n.(n.+ l)*A(n.+ n.,f.+ f.+ 2,m.+ ra.+ 2 ) 

J J 1 J 1 J 1 J 
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+ 



1) •A(n^ + 


n .+ 
3 


2,X .+ 
1 


X . ,rr 

3 


1 . + m . + 
1 3 


• 2 ) 


A(n.+ 


n .+ 
3 


2 ,X .+ 
’ 1 


X .+ 
3 


2 ,ra^ + 


m .+ 1 ) 
3 


• B( n^ + 


n .+ 
3 


2,J^ .4- 


X .+ 
3 


3 + 


m. ) 


. • B(n . + 
J 1 


n .+ 
3 


3,^ 


X .+ 
3 


2 , m . 4- 

1 


"3 ^ 


m .n . ■ B( n . + 
3 J 1 


n .+ 
3 


1 , + 
1 


X .+ 
3 


3,m^+ 


"3 ^ 


m .1 . • B(n. + 
D D 1 


n .+ 
3 


3,^ .+ 


X .+ 
3 


1 ,m . 4- 
1 


m. ) 


'A(£^ + 


i .+ 
3 


2 , n . 4- 
X 


n .+ 
3 


2 , m . 4- 

’ 1 


m .+ 2) 
3 


[J 20 ;(jij+ l)-z3 •A('^£ + 


i .+ 
3 


1 ,n . 4- 
1 


n .+ 
3 


2 ,ra . 4- 

1 


m . + 2 ) 
3 


L'^ia(n^+ l)-z] *A(X^ + 


i .+ 
3 


2 , n . 4- 
1 


n .+ 
3 


1 ,m . 4- 

1 


m . + 2) 
3 


. • A(^ . + 
J 1 


i .+ 
3 


2 , n . 4- 


n .+ 
3 


3 5m^4- 


m. ) 


. -A(.£ . + 
J 3 


Si .+ 
3 


3,n^ + 


n .+ 
3 


2 ,m . 4- 

1 


) 


m.(m.+ 1 .+ n .+ l)*A(i.+ 
3 3 3 3 ' ' 1 


1 .+ 
3 


2 ,n . + 
1 


n .+ 
3 


2 , m . 4- 

1 


m. ) 


.(i .+ 1) -A(je .+ 


X . , n . + n . • 
J 1 J 


t-2,m.+ m.+ 2 ) 

13 


'inj(rij+ 1) -A(£^ + 


X .+ 
J 


2 , n . 
1 


n . ,m .+ m .+ 2 ) 

313 


A(^.+ 


X .+ 
J 


2 ,n . 4 - 
1 


n . 4 - 
J 


2 , m . 4- 

’ 1 


m .+ 1 ) 
3 


• E{i ^ + 


X .+ 
J 


2 ,n . 4- 
1 


n .4- 

J 


3,m^+ 


™3 ^ 


• B(i. ^ + 


X .+ 
3 


3,n^+ 


n .4- 
J 


2 , m . 4- 

1 


"3 ^ 


ra .£ . ' B(i> . + 
J J 1 


X .+ 
3 


1 , n . 4- 
1 


n .+ 
J 


3,m^ + 


"j ^ 


m .n . • B(J^ . + 
3 3 1 


X .+ 
3 


3,n^ + 


n . 4 - 
J 


1 ,m . 4- 
1 


m. ) 


2 

+ *A(n^ + 


X .+ 
3 


2 ,X .+ 


n .4- 

J 


2 , m . 4- 

1 


m .+ 2) 
3 


l)-z] *A(n^ + 


X .+ 
3 




n .+ 
J 


2 , m . + 
1 


m .+ 2) 
3 


[i^(n^+ 1 ) -z3 -A(n^ + 


X .+ 
3 


2 ,X .+ 
1 


n .4* 

J 


1 , m . 4* 

1 


m .+ 2) 
3 
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+ 


^ra . *A( n . + 


i .+ 


2 ,£ .+ 


n .+ 3 J ni . + 


ra . 


) 




J 1 


3 


1 


J 1 


3 




+ 


. • A ( n . + 


£ .+ 


3,^ .+ 


n . + 2 , m . + 


ra . 


) 




D 1 


3 


1 


3 1 


3 




_ 


m.(ra.+ n .+ £ .+ l)*A(n.+ 


£ .+ 


2,£ .+ 


n .+ 2 , m . + 


ra . 


) 




3 3 3 3 ' ' 1 


3 


1 


J 1 


3 




_ 


h^.(£.+ 1) -A(n + 


£ . ,£ 


. + n 


- 2 , m . + m 


* 2 


) 






3 


1 J 


1 3 








hn.(n + 1) -A(n.+ 


^ . + 2 


!,i.+ n.,m.+ m.+ 


2 


) 




3 3 ' ' 1 


3 


X 


3^3 






+ 


A(n^ + 


£ .+ 
3 


2,£ .+ 
1 


n .+ 2 , m . + 

J 1 


ra .+ 
3 


1) 




. • B(n . + 


£ .+ 


2,£ .+ 


n .+ 3 ,ni . + 


ra . 


) 




3 1 


3 


1 


3 ^ 


3 




_ 


. • B( n . + 


.+ 


3,£.+ 


n .+ 2 . + 


ra . 


) 


* 


3 1 


3 


1 


J 3. 


J 






ra .£ . • B( n . + 


^ .+ 


1 ,£ .+ 


n .+ 3 J ^ 


ra , 


) 




3 3 1 


3 


1 


3 1 


3 






ra .n . ‘ B( n . + 


i .+ 


3 


n .+ 1 , m . + 


ra , 


) 




3 3 1 


3 


1 


3 1 


3 




- 


‘A(£.+ 


n .+ 
3 


2 , n . + 
1 


. + 2 , m . + 
J 1 


m .+ 
3 


2) 


+ 


[ijy(n^+ l)-z] 'A(^j + 


n .+ 
3 


l,n.. 


i« .+ 2 , m . + 
J 1 


m .+ 
3 


2) 


+ 


[ija(£ .+ 1 ) -z] ‘A(£ .+ 


n .+ 


2 , n . + 


X . + 1 , m . 


ra . + 


2) 




~J 1 


3 


1 


J 1 


3 




+ 


?sam . *A(£ . + 


n .+ 


2 ,n + 


^ + 3 , ra + 


ra . 


) 




3 1 


3 


1 


3 ^ 


3 




+ 


i-^ra •A(£ + 


n .+ 


3,n + 


j 6 2 , m . + 


ra . 


) 




3 1 


3 


1 


3 1 


3 




_ 


m.(m.+ n .+ £ .+ l)'A(jf..+ 


n .+ 


2 , n . + 


J? . + 2 , ra . + 


m . 


) 




3 3 3 3 ' ' 1 


3 


1 


J 1 


J 




_ 


(n + 1) *A(£ . + 


n . , n . + i 


2 , m . + m .• 


■i- 2 


) 




3 3 3 


J 


^ J 


1 3 






_ 


l)'A(i.+ 


n .+ 


2 , n . + 


. ,m . + m .• 


+• 2 


) 




3 3 ' ' 1 


J 


1 


J 1 J 








• A(^.+ 


n .+ 
D 


2 , n . + 
1 


. + 2 , ra . + 
3 1 


ra .+ 
3 


1) 




Jgam . -Bfi .+ 


n .+ 


2 , n . + 


£ .+ 3,ra.+ 


m . 


) 




3 1 


D 


1 


3 1 


3 




_ 


^^ra . • B(X . + 


n .+ 


3 , n . + 


. + 2 , m . + 


ra . 


) 




3 1 


D 


1 


J 1 


3 




+ 


ra .n . • B(^ . + 


n .+ 


1 , n . + 


.+ 3 J ra . + 


ra . 


) 




3 3 1 


D 


1 


J 1 


3 




+ 


ra .i . • B(£ . + 
3 3 1 


n .+ 
J 


3,n^ + 


i . + 1 , m . 
J 1 


ra . 
3 


)} 
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and 




A(n.+ n.+ 2,^.+ £> ,+ 2,m.+ ra.+ 2) 

' 1 j 1 j 1 j ' 

A(jC.+ i'.+ 2,n.+ n.+ 2,m.+ m.+ 2) 

1 J ’ 1 J 1 J ^ 

+ { A(n.+ ^.+ 2 ,jC.+ n.+ 2,m.+ m.+ 2) 

— J 1 J 1 J ' 

A(jC.+ n.+ 2,n.+ jC.+ 2.m.+ m.+ 2) } (A-5) 
13 ij ’13 ' ' ^ 



Note that if the expression (54) is used to evaluate the A’s 

and B*s that each of the terms in (A~4) and (A-5) above must be 

. . V +X +u» 

divided by Of 

2 • <Sl (H - E) and <C| (H - E)|\ > 

Each of these matrix elements is made up of 26 terms which in 
turn are made up of an integral of the form (53 a-d) where A*s and 
B’s are as defined in ( 44 ) to ( 47 ) and immediately following. 

Now letting \ in (A-3) above be char act or i;^ed by n,f,m, <S ,C | ( H-E) |\ ^ 
is 



<S,c| (H-E)lx > 



= - E) 'A 

+ [? 2 a{i+ 1) -z] -A 
+ [^ga(n+ 1) -z] -A 
+ -A 

+ ?§0;ro • A 

-m(ra+ i + n+ 1 ) • A 
1) -A 
%n ( n+ 1 ) • A 
+ A 



(£ + 
s , c ' 


1 ,n + 


2 , m + 


2) 


s,c^^ 2 ,m+ 2 


) 


s ,c' 


1 5 n + 


1 , m+ 


2) 


(£ + 
s ,c' 


1 , n + 


3 ,m 


) 


(£ + 
s ,c' 


2 , n + 


2 ,m 


) 


(£ + 
s ,c' 


1 , n + 


2 , m 


) 


s ,c' 


1 J n + 


2 5 m+ 


2) 


(£ + 
s , c ' 


1 , n , m+ 2 


) 


(i' + 
s , c ' 


1 , n + 


2 , m + 


1) 
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l,n+ 3,m ) 



(J^+ 2,n+ 2,m ) 

s « c 



J0m« B ('^ ,n+ 3 >ni 
s , c 



) 



+ mn*B (ji+2,n+l,ra ) 

s j o 

+ E) *A^ 2,m+ 2) 

+ C^(n+ l)-z]*A (n,-^+ 2 ,m+ 2 ) 

S j o 

+ [igX(Ji+ 1)-Z]*A (n+ 1,£+ l,m+ 2) 



s , c 



^ra-As ^(n+ 1,X+ 3,m ) 



^^Oon*A^ 2,J?.+ 2,ra ) 



-ra(ra+ £+ n-** 1)-A^ 2,m ) 

^n(n+ 1)*A^ 1 ,'^'^ 2,m+ 2) 



-2^(j£+ 1)-A^ l,^,m+ 2 ) 



s , c 



( n+ 1 + 2 , m+ 1 ) 



^ra*B^ ^(n+ 1,£+ 3,ni ) 



-/Vra-B (n+ 2 ,£+ 2,ra ) 
s j o 



mn • B (rj ,jC+ 3 jHi 
s , c ' 



inX*Bs 2,£+ l,m ) } 



(A-6) 



where it is assumed A , B are chosen when the left hand side is 

s s 

<S|(H-E)|x> and A^ , B^ when it is <c|(H-E)lx>. 

3 . Determination of the Required Integ rals 

Clearly not every term in the above matrix elements requires a 
different integral. Many of them are identical or have ;?cro 
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coefficient. For example, for the 34 element basis set (M = 5 ) , 
while there are 595 independent elements each in H ^nd S, repre- 
senting 52 terras and 4 terms each, respectively, or a total of 
more than 33,000 terras, only about 700, or slightly raore than two 
percent, use different integrals and all 700 different integrals 
can be generated using the recursion relations proved in Appendix B 
by actually carrying out the evaluation of about 100 integrals for 
which /i = 2 . Similarly, in evaluating <s| (H-E)lx^^nd <c|h-H)|x>, 

where there are 34 elements in each expression, each element in turn 
$ 

made up of 26 terras, for a total of 1768 integrals required, only 
268, or slightly less than one sixth are distinct and again, using 
the recursion relations only IO3 terms for jl - 2 need be evaluated 
explicitly and they can all be found from 36 different numerical 
integrals. Nevertheless it is the numerical evaluation of these 
integrals which overwhelmingly dominates the time required to solve 
the complete problem and hence the requirement only to evaluate 
those terms needed. 

To insure that no terms were missed in computing the integral 
tables, preliminary computer programs were written to examine each 
of the matrix elements and to display in tabular form which integrals 
were needed to construct the matrix elements for a given basis set 
size. Tabulated were all distinct integrals required for which the 
coefficient was non-zero in at least one case plus those required 
in order that the recursion relations could be used to generate 
the terras required for /i = 2 . From this information the integral 
table programs were constructed so that only those terms required 
were actually calculated, and none that would be required was omitted. 
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APPENDIX B 



Integral Recursion Relations 

1 . R ecursion Relations for [i ^ 2[ 5l] 
Consider the integrals 



A(v,X,Li)= yr^dr^r^ 



(B-1) 



and 



B(vA,u)= 



(B-2) 



and rewrite (B-l) in the form 



A(vA,u)= 4 - 2'i>^2'=°='‘i2)''("i''2>*1*: 



using the law of cosines. This transforms immediately to 



A(V,X,U')= A(V+2 ,X ,M'-2) + A(V,X+2,U-2)- 2B( V + 1 ,X +1 .M' . 



(B-3) 



Similarly B(v,X>M') becomes 



B(v,XjU)= B(v+2 ,X ,b-2) + B(v,X+2,a-2)- 2^r r ^ ~'^F ( r ^ , r ^ ) cosQ^^dr ^dr . 

Let the integral in the third term be designated I, then 



I = A(v+l,X+l,U'-2)-^r^"^r2"^r^2^F(r^,r2)si 



sin «i2^"i'^"2' 



Now recall that the angular part of the dr^dr^ integration is 
simply sin0^2^^]L2’ angular part if the second terra in I is 






^ ~4 

(l‘- xl- 2r^r2COsej2)’^ ®i"^®12‘’*12 



= _)_ _i_ cr 

f-2 jL 



(ij+ r^- 2rjr2Cosej2)"2' jsin-9^^09^^ 



\x-2 



1 . 2, 



d9 



12 
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which raay be integrated once by parts to give 






H-2 yi2^ 



and therefore 



I = A(V+1 ,X+1 ,U-2)+ B{vA,t^) 

1-L 



from which it follows immediately that 

B(v,X.l^)= {b(V +2 ,X ,U'- 2 )+B(v,X +2 ,ii- 2 ) - 2 A(V + 1 ,X+1 ( B-4) 



2. Recursion Relations for u = 1. 



When u = 1, (B-1) becomes 



A(V,X,1)= dr, dr 



12 



1 2 



Now recalling that in the coordinate system in use ^ ^ 2 ’ 
l/r ^2 can be expanded in a series of Legendre Polynomials 



= T. 



P. (cos9_) 



"l2 ^=0 



and the angular integration becomes 

I 






CO r 

a 



But P^(u) = 1 so the integral on the right can be rewritten 



r.i 



and by virture of the orthogonality of the Legendre polynomials, 
all terms in this expansion vanish except for f = 0. Therefore 



A(V,X ,1) = A(V-1 ,V,2) . 



(B-5) 
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In a similar exj)ansion in the case of B(V,X>1) terms but 

f = 1 vanish, giving 



B(V,X,1) = -| A(V-2,X+1,2) . (B-6) 

3 . Recursion Relations for A (v,X,M/) and B (v,X,ll) 

s , c s , c 

Referring to (53a) 

Ag(V,X,VO = A^(V,X,M.)+ A^CV.X.M-) 

and applying (B-3) above 

A^(V,X>U')= A^(X+2,V,U-2)+ A^(X ,V-2,p.-2)- 2B^ (X +1 , V + 1 ,M. -2 ) 

+ A2(v,X+2,a-2)+ A^(v+2 ,X ,p.-2) - 2B^ ( V+1 ,X +1 jLi -2 ) 

= Ag(V ,X+2 ,M.-2) + Ag(V+2,X >^‘’-2) - 2Bg(V+l ,X +1 )M'-2) 



and similarly with B (v,X)M>). Thus A and B satisfy the 

s s , c s , c 

normal recursion relationships for [i > 2 . For \x = 1 

Ag(v,X>l) = A^(X,V,1)+ A^CvA.I) 

= A^(X-1,V,2)+ A^Cv-l.X-S) 



but 



so 



A2(V-1,X,2) = A^(v-l,X,2)« A2(X,V-1,2) 



A 



V ,X 5 1 ) A^ ^(V“1,X,2)+ A^ 2(X”1>V,2) - A^ 2(XjV~1,2) ( B- / ) 



where A. is chosen when A is to be evaluated and A^ is chosen when 
is 3 

A^ is to be evaluated. For B^ ^(v,X,l) the same procedure yields 
^(v,X,l)— “ 3^^ 



(B-8) 



which completes the derivation of all the relations needed to 
evaluate the matrix elements listed in Appendix A. 
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APPENDIX C 



Ray leigh -Ri tz Variational Method 

Let an arbitrary quadr atically integrable wave function Y be 
expanded in terms of the eigenfunctions of a bound system 



whet e 



T = ? a.Xj 



(H-E)IXi> = 0 



then on the function Y the expectation value of H is 



<Y I H| T> = E a . a ,<X . n| X > 
i,j 1 J 1 ' J 



= E a.rE. . 

i " " 



If E is the ground state of the system, then E ^ E. so 
o o 1 



E E|a.|^ S E la.I^E. 

Ol ' 1 ' H 1 1 



whence 



< V I H I Y > s E 



(C-1) 



(C-2) 



assuming 



But 



< Y Y > = 1 . 



< y1h|y > = Ey 

the expectation of H on Y , so if Y is found so as to minimize E^ 
then E^, > and E^ provides an upper bound on the value of E^ . 
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In general 



= < yIhIy >/<y1'1' > . 



(C-3) 



If now a trial u^ve function is chosen and expanded on a suitable 



basis 71 . 

' 1 



N 



Y = E b.n . Y 

t 1 ' 1 



(c-4) 



whence 



S b% <?7 I h|? 7 . > S b?b.H.. 

. . 1 J 1 3 . . 1 J IT 

= 



t * I 

E b . b .<71 . 11 . 

i,j 1 J 'i' '3 



E b.b .S. . 
i,j 13 iJ 



where H. . = < 71 . | h|t? > and S. . = < 1? . jl? .> . 

ij ' X j xj ' X j 

to satisfy 



Now choosing the b's 



bb. 

1 



= 0, 



i = 1, . . . ,N 



gives 



bE 



t 



Ob. 



1 



s 

1,3 



-K- 



b.b .S. .+ E 
1 3 13 



t 




b% .S . . 
1 3 13 



b 

bb. 



S 

i.3 



b*b .H. . 
1 3 13 



from which arise the N linear equations 



E(H. E^S. .)b. = 0, 

i' xj t xj' X 



3 = 1 



,N 



(c-5) 



which is just the first part of the Harris expansion technique, 
wherein the smallest eigenvalue gives an upper bound on 

The above method extends readily to the excited states [ 5^1 • 

This is achieved by choosing a trial wave function which is 
orthogonal to all the states lying below the desired one. Such a 
function must take the form 

N-1 

Y = cp - E X„<X^l9 > (C-6) 

n=0 n 
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where cp is arbitrary and the are the lower lying eigenstates. 

Then expanding Y according to (C-1) it is seen that the a^ are 

zero for i less than N, the index of the state desired, from which 

it follows that ^ . 

t N 

Hylleraas and Undheim [32] have shov/n that this condition is 
automatically met by (C-5), the n^^^ eigenvalue of that equation 
being greater than or equal to the energy of the n^^^ level of the 
system. Thus it is seen that an important bonus in the Harris 
expansion, method is a cataloging of the energy levels of the bound 
system formed by the scattering atom and the incident electron. 
Included in this cataloging are the autoionizing levels of the 
bound system which give rise to many of the observed scattering 
resonances [ 44 J. 
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APPENDIX D 



Discussion of Numerical Calculations 



1 . The Eigenvalue Problem 

From Section II C.l, a solution for eigenvectors and eigen- 
values is sought for the equation 

(H - XS)c = 0 (D-1) 



Since S is not the identity matrix (D-1) is not in the form usually 
associated with eigenvalue problems occurring in physics where the 
space is orthonormal and S reduces to the identity matrix. This 
generalized eigenvalue problem has been extensively discussed in 
the mathematical literature [4l]- In the present case, however, 
it is possible to reduce the problem to standard form and treat 
it by the simpler methods available. 

Consider first the matrix S. Its elements are the inner pro- 
ducts of the chosen basis vectors 



S, ^ = <X J\' > . 






1 J 



(D-2) 



Since the basis functions are real S. . = S.. and therefore there 
exists a unitary transformation which diagonalizes S CSSl; Physi- 
cally this is equivalent to transforming to a nev/, orthogonal set 
of basis vectors. In this new set 







<X'.|X!> 
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and 0 for all i. But diagonalizing a matrix in this manner 

means that the new matrix has as its diagonal elements the eigen- 
values of Therefore it can be concluded that all the eigen- 

values of S are positive and hence S_ is positive definite. 

Now consider the positive definite quadratic form associated 
with X -^S'X. There exists a real non-singular transformation 
^ and a vector y such that x = U and such that (Ref. 53, p. 337) 



y 




• y 



- 

y -y 



called the '"cannonical form’ 



of X -S *x , 



where U 



(D“3) 

(U Hence 



U‘^.S-U“^ = X (^’"4) 

where X identity matrix. So 

S = U^U . (D-5) 

In other words, a positive definite symmetric matrix can always be 
factored into the product of a real non-singular matrix and its 
transpose. But because S is symmetric, this factorization is not 
unique, for 

N 

S . . = 5: H, . U. . 

k=l kj 



represents n equations, only n(n+l)/2 of which are independent. 

Therefore one is free to choose n(n-l)/2 of the U/, arbitrarily, 

' ' Urn 



3 

Unless one of the nev; basis vectors is null, which would mean 
the dimensionality of the space is reduced - a possibility which 
can be excluded from the present problem. 



89 



and it is convenient to choose 



U. =0 £ > m 

xra 

for then U becomes an upper triangular me\tr ix (one in which all 
elements below the principal diagonal are zero). 

With the above choice of factorization of S^, (D-1) can be 

r eivT: itten 

Xu"^*u)u"^ -U'C = 0 



and the original problem reduces to 



- X^)-x = 0 



(D-6) 



- - . . ..... -T - 1 

where x = U-c. Since H is symmetric it is trivial to show U -H-U 



is also symmetric, for 



. = E uT^h. .u 7.^ - s u7^h, .u7 

'mj ik IJ ik im ki j 



-T 

ik 



= E u7^h .u"^ = (u"^-iru"^) . 

i ,k Jk kj Im '' jm 

and the problem has been reduced to the standard problem of a real 
symmetric matrix in an orthonorraal space. 

The advantage of the form of the factorization of S chosen is 
that because U is triangular its inverse is particularly easy to 
obtain. The factorization is done by the square root method of 
Cholesky C54l and U ^ is then found by solving the algebraic 
equations [55] 



u7^u. . = 1 

11 11 



i = 1, , . . ,N 



E u7 .U., = 0 , i< k 
j=l 3^ 



(D-7) 



for the U.^. . 
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Because the range of values of the elements of S and H v;as 
often several orders of magnitude (sometimes as many as fifteen 
or more) an additional check on the accuracy of U ^ as found by 
the above method v.^s used. If the U ^ found above is taken as a 
first approximation to the exact inverse, then the following 
iteration scheme can be used to improve the value of the inverse L 56 ] . 

Let be the i approximation to the exact inverse, U , and 
define the error matrix 

' R. = I - C. -u . (D-8) 



Then 



and 



wher e 






C.+ R.C. = C.+ (I- C. • U) *C. 
— 1 — 1—1 —1 — —1 — ' — 1 



R . , 
—1 + 1 



i ■ (L ■ ‘U 



(1 - c. -U) = 



r" 

—I 






R = I - C *U 
— o -- — o — 



(D-9) 



and C is taken to be the value of U ^ found by the Cholesky 
method above. Clearly, if the norm of R is less than one the 
method converges quite rapidly, and since the value of is 
already quite accurate, i • In practice no more than one 

iteration of C was ever required to make the elements of R .< 10 
The eigenvalue problem was solved using the Jacobi variable 
threshold method [4l,55l* This method proved to give the most 
accurate results although it is known to be considerably slower 
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than other methods tried. The accuracy of the eigenvalues was 



verified by the following method. For a given eigenvalue, X , the 
determinant |H“X^s 1 was formed and compared with the determinant 

“7 

when X^ was replaced by ^ where 6 was typically 10 (compared 

to X which was on the order of one). Then X was considered cor- 

M' 

rect to six places if the respective determinants were related by 



H-(X 1 S)s|<iH - X^s|<|h - X + 6)S| 



(D-10) 



and the first and last terns above were of opposite sign. Accuracy 
of the eigenvectors was not so clearly established, though a 
similar procedure to the above was followed. The product (H-Xs)c 
was formed and it was noted that as X was perturbed as above, all 
elements in the resultant vector changed sign. IVhile this indicated 
that the eigenvectors probably are reasonably accurate no quanti- 
tative bound could be placed on their' error. Further indication of 
the accuracy of the eigenvectors can be inferred from the close 
agreement of the phase shift results for hydrogen with the accepted 
values found by Schwartz [l2]. 

2 . Numer ical Integration 

From the beginning it was felt that the simplest and most 
reliably accurate means of evaluating the necessary integrals 
numerically would be some form of Gaussian quadrature C57J- Any 
numerical integration scheme attempts to replace the integration 
by a finite sum so chosen as to make the error tolerably small 




w . f ( X . ) 



(D-ll) 
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where the w^ are weights to be assigned to the function evaluated 

at the sampling points x^, in other words a weighted average of the 

function at the various sample points. If the are fixed, equally 

spaced points then the n + l can be chosen to define a polynomial 

of order n with which to approximate f(x), and if f(x) is itself a 

polynomial of order n or less, the integration will be exact. If 

however, one is free to vary the x^ also, there will now be 2n + 2 

parameters available and it will be possible to define a polynomial 

of degree 2n + 1 with which to approximate f(x). Thus, given the 

values of x. and w. in each case, the same amount of labor results 
1 1 

in a more accurate integration in the second case. 

Adopting this approach, one can approximate f(x) with a 
Lagrangian interpolating polynomial of degree (n + 1) of the form 



f(x) = g(x)F(x)= E L (x)g(x)F(x )+ R (x) 

1=0 1 in 



(D-12) 



v;her e 



L.(x) =7T : 



(x-x^.) 



(D-13) 



j=0 (x . -X .) 

■a- 1 3 

and R (x) is the remainder term, given by 



-TX (5) g(x)/(n+l) : , a<5>b (D-l4) 



i=0 



The reason for factoring f(x) into g(x)F(x) will become apparent 
below. Now if F(x) is a polynomial of degree 2n + 1 then F^^ ^^5) 

must be a polynomial of degree n. Letting 



p(^ 

(n+l ) 



ill 
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where ^ polynomial of degree n, (D-l4) becomes 



- <1A^)TT 

i=0 



Integrating (D-12) gives 



Now (D-16) is of the form (D-11) where 



(D-15) 



pb n pb pb n 

\ g(x)F(x)dx= E ^F(x )\ g(x)L (x)dx + \ g(x)q (x);i (x-x,)dx ( D-l6 ) 

" i=0 ^ 



p u n ^ X - X . 

^ j=0 

jVi 



n (x-x^.) 

r^'T-x ;■)' 

1 J 



dx 



(D-17) 



except that the x^ have not been chosen. The object is to so 

choose the x. that the second term of (D-l6) vanishes. This is 

n 

easily done if q (x) and J~f (x - x.) are e<panded in a series of 
" i=0 ^ 

polynomials orthogonal on [a ,b] with respect to the weight 
function g(x). A weight function appropriate to a particular 
set of orthogonal polynomials is often a natural factor in the 
function to be integrated and therefore the particular expansion 
to be chosen may be dictated by the problem, hence the utility in 
carrying the factor g(x) through the derivation. 

Let the polynomials chosen for the expansion be designated 
P^(x) whereupon 

(D-18) 



n n + 1 

^ ‘’‘■’‘d ° !-0 

i=0 



and 



q (x) = E c P (x) 
1 m=0 tn 



(D-19) 
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Combining these two expansions in (D-15) and applying the ortho- 
gonality condition 

.b 

\ g(x)P^ (x)P^(x) dx = 0, i / m (D-20) 

a 

reduces the error terra to 




g(^) %(-'<) 



Tt 



i=0 



n pb 

(x-x.)dx = E b.c. \ g(x)[p.(x)3 dx 

j=0 ^ J ^ 



(D-21) 



Now if the x. are chosen to be the zeroes of P ,(x), then the 
1 n+1 ' ' 

expansion (D-l8) reduces to 

7T (x-x. ) = b P (x) 

^ 3/ n + l n + 1^ ^ 

i=0 

i.e.j b^. = 0 O^j^n, whence (D-21) vanishes, and the integration 
is exact. Thus evaluating the function at n+1 carefully chosen 
points suffices to integrate it exactly if it can be exp^ressed 
in terms of a polynomial of degree 2n + 1 or less, a considerable 
improvement over the equally spaced interval methods, as was 
promised above. The one remaining problem is, then, to evaluate 
the wjs and find the x.’s. Fortunately, this has been done for a 
large class of orthogonal polynomials by Stroud and Secrest [58]. 

Although the integration limits and the form of the integrands 
( 49) j ( 5 l) > ( 56 ) and (57) suggest the use of the Laguerre poly- 

nomials for the numerical integrations, the almost periodic nature 
of the functions and the ever increasing amplitude in the absence 

— X 

of the e factor dictated the use of a finite interval formula 
where several intervals could be integrated separately and then 

— X 

added together to give the final result. Since the factor e 

forces the integral to zero at large x, examination of the integrands 
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indicated that an upper limit on the integration of x = 100 would 
introduce a negligible error from the neglected region in the worst 
possible case. Because the variable change to conform to the con- 
vergence limits of the polynomials used in the expansions (D-l8,19) 
is simple and because the weight function is unity the decision 
was made to use Legendre polynomials in the integration. 

Because of the almost periodic nature of the integrands, when 
the instantaneous frequency of 1he sinusoidal variation is large the 
areas under the integrand above and below the axis are nearly equal 
and substantial loss of significance occurs when these two values 
are subtracted. Often as many as eight to twelve significant 
figures can be lost due to subtraction, hence great care must be 
taken to avoid the problem. The approach adopted was to seek a 
function whose exact integral is known and which approximates the 
desired integral closely enough so that the difference between them 
is small, then evaluating this small difference numerically and 
adding the result to the known value of the approximating function. 
This method can be compared to the measurement of the difference 
between two large quantities, say frequencies. Often the difference 
can be determined quite accurately even though the absolute magni- 
tude is known only approximately. 

In executing this approach to the integration of the functions 
at hand the first step is to find all the zeroes of the integrand 
between arbitrarily chosen minimum and maximum values of x. These 
values are then used as the end points of the subintervals in an 
integration of the form (D-22) 

n I ^ rr 

I(^)=\ f(x)dx+ ?_ |\ 

i 



^i+1 

l^f ( x) -h^ ( X h^(x)dxj'+^ f(x)dx 
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where the are the n + 1 zeroes of f(x) found above and the h.(x) 
are functions whose exact integral is known and which approximate 
f(x) over the interval indicated. The h^(x) are chosen as indicated 
in Figure 23 . Note that h^(x) is amplitude modulated from segment 
to segment so that f(x) - h^(x) is always mostly positive and 
smaller than f(x) by at least an order of magnitude regardless of 
the sign of f{x). It is this difference function which is inte- 
grated numerically, and since h^(x) is a different function for 
each half cycle of f(x) the numerical integration is done in 
segments between zeroes of the function using a 32-point Gauss- 
Legendre quadrature on each segment. 

The form cliosen for h.(x) was 

1 ' 

h^(x) = (1 + “ )sin(A^+B)e ^x^(l-e (D-23) 



where p is +1 if f(x) is negative and -1 if f(x) is p^ositive, 

B = 1 if I or I is to be found and B = 0:/2(Z'Ky) if I or I is 

^ sc ^^^13 

to be found, and 



A = —V . (2 * 



(D-2Zf) 



i + 1 i 

SO chosen that sin(Ax + B) vanishes at the end points of each 

segment and nowhere in between. 

When chosen in this manner the integral of h^(x) alone is 

- Px 3 

exact and has a rather simple form. If (1 - e ) is expanded 



four integrals of the form 

-a^x 

~ J x^clx, X = 1,2, 3,^4) 



(D-25) 



Z . 

1 
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Figure 23 . 

Example of Metliod of Choosing h^(x) 






result, where is 1 , 1 + p , 1 + 2p and 1 + 3p in turn. These 
integrate immediately to 



s, = y 



I Y, 5 in ( jO ) / 

D=i (-o)J(y-^i-j) : ^ 



_"^f^i+l y+i-j 
^i+l 









(D- 26 ) 



where 0 and p satisfy the relationship 



+ (-l)^A = p(cos© + (-l)^sin9) 



and finally 






i+l 



z . 
1 



h^ ( x) dx 



10^ ^^l”^^2^ ^^3~ ^4^ 



(D-27) 
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Expr ess ions ( D-26 , 27) arc valid for 0. During the course of 
evaluating the phase shifts, however , it is necessary to evaluate 
the various integrals for y = -1. For this case a slightly dif- 
ferent scheme must be used to evaluate h^(x). If one replaces 
(1 - e ^^)/x in (D-23) by the integral identity 




then 



(D-28) 



h^(x) = 



P(1 



" So""" 






dy 



(D-29) 



Carrying out the x integration first and making the appropriate 
variable changes yields 



p(f +1)+1 



- z . w - z . ^ w 
O " t o . 

2 5 dw 

w + A 



(D-30) 



and finally 

r ^ p 

\ h.(x)dx = “ P(1 + — )-A‘(S -2S +S^) 

J ^ 10^ ^ o 1 2^ (D- 31 ) 

i 

Unfortunately (D-30) must be integrated numerically also, however 
it is a relatively slowly varying function of w and is always 
positive so integration by a 32-point Gauss -Legendr e quadrature 
gives more than adequate accuracy. 

To complete the integration of (D-22) the first and last terms 
were evaluated using single 512-point Gauss -Legendr e quadratures. 

Since the accuracy of the numerical integration is crucial to 
the accuracy of the phase shifts calculated for helium ion scat- 
tering, a great deal of effort was expended to insure that the 
methods were both correct and internally accurate. 
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The first and most obvious check is to increase the number of 



sample points in the numerical integrations. This was done by 
increasing the quadrature method on each interval to 64 points. 

t h 

Changes to the integrals on doing this were typically in the 12 
t h 

to l4 digit, indicating that the differences were principally 
due to round-off error since the maximum number of digits available 
was l6 . 

The integration scheme described above grew out of a similar 
method which was used extensively prior to the present one for 

I 

the same calculations. The present method was adopted in an effort 
to extend the useable range of the parameters involved in the 
integrands. The problems with the former method arose principally 
because it was not sufficiently accurate when the oscillation fre- 
quency of the sinusoidal function was high. However in the region 
where both schemes wore adequate (which proved to bo almost every- 
where) the integrals were consistently r epr oduceable to 12 to l4 
places. The choice of integration intervals and of the h^(x) in 
the earlier method was sufficently different that it constitutes 
an almost independent chock on the internal accuracy of the 
integrals . 

Three different methods were used to make the integrals exactly 
intograble so that the numerical results could bo compared directly 
v/ith known answers and to provide a check on the correctness of the 
expressions. The first was to let k ^ while holding the ratio 
k/oc constant. This forces the logarithmic terms to zero and re- 
duces the integral to a standard form which can be integrated 
directly. The direct evaluation was done on a Hewlett-Packard 
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Model 9 SIOA computer which gave results to ten significant figures. 
These results were then compared with the numerical results and 
except in the cases where k/o. was at its extreme limits, agreement 
was obtained to eight to ten places. This method was tried only 
with the older integration scheme and so some decay at the limits 
of k/oc was expected. In addition, the transition from the form 
of the integral when k was finite to its form as k appeared 

smooth and with no obvious discontinuities or other strange behavior 
in the transition region. 

A second check took advantage of the fact the computer program 
was written to allow for variation of the nuclear charge, Z. 
Although the intent was to limit the value of Z to integers, 
nothing in the program or the mathematics forbade varying the 
charge continuously from one to two. When this was done, it 
was again seen that the result for Z = 2 transitioned smoothly 
into that for Z = 1 with no discontinuities as the charge was 
var ied . 

Since for Z = 1 the integrals all become exactly integrable, 
the computer program did not do the calculations required for the 
numerical integration in the case of Z = 1 , but skipped to a sub- 
routine which evaluated directly the exact integrals. The final 
comparison method was to circumvent the exact integration routine 
for Z = 1 and do these integrals numerically, using the method 
described at length above. Again the agreement between the two 
calculations was excellent, the integrals and the phase shifts 
typically being identical in the first twelve to fifteen places. 
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J 



The final check was to lest the sensitivity of the entire 
calculation to perturbations in k. It was noted that even at the 
extreme limits of k/o: where accuracy was expected to be poorest, 
only some of the integrals showed evidence of loss of significance 
from the subtraction problem mentioned above. For these integrals, 
however, a small perturbation in the value of k may have a rela- 
tively large effect on the integration and if the inaccurate 
integrals dominated the computation then they could significantly 
distort the results. As before when this perturbation was arti- 
ficially fed into the problem its effect was well within the limits 
considered acceptable (results unchanged to six places or more for 
a perturbation of k in the seventh place) except at the extremes. 

Since none of these tests indicate any serious problems with 
accuracy it is felt that the errors inherent in numerical inte- 
gration have been kept well within tolerable bounds in this problem 
and the results can be relied upon with considerable confidence. 

It also proved to be true that those regions where accuracy was 
lost were easily avoidable. 

3 . Evaluation of r(l - ia) 

From the definition of the gamma function 




CO 



(D-32) 



0 



one finds 




CD 



r(l-ia) = \ G X -^“dx 
but e X, so (D-33) can be rewritten 



dx 



(D-33) 





= A + iB 
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whence 



atg [r(i-ia)} = tan ^ — 

A 



(D-35) 



and the evaluation of the complex gamma function is reduced to 
evaluating two real integrals. The integrals have to be evaluated 
numerically, however, and are not suitable for this in the form 
(D-34) because the rapid oscillations of the integrand near the 
origin could contribute to a substantial error in the numerical 
integration. To overcome this problem note that an integration 
by parts of each of the integrals (D-34) introduces a factor x 
into the integrand and this has the desirable effect of reducing 
the contribution of the integrand near the origin. Each re- 
petition of the partial integration introduces an additional 
factor of X into the integrand. Carrying out this process three 
times transforms the integrals (D-34), after considerable algebra, 
into 



B = 

and 
A = 



-a 



(l+a^) (Zf+a^)(9+a^)‘^0 



-a 

( l + a ^)( 4 + a ^) ( 9 + a ‘^)'^0 



CD f 

^ x^e ^ "^^ ( a^ -1 ) s in (a''":x ) - (a"^-ll ) cos ( a-''*>c) j*dx 

(D- 36 ) 

2 ^ x^e ^ ‘^(a^-11 ) s in (a'^T^x) + ^(a^-l)cos(a.‘'":x)j‘dx 

(D-37) 



which forms are now quite suitable for numerical integration. 

The integrations of these functions were carried out using a 
256 -point Gauss -Legendr e quadrature over 15 arbitrarily chosen 
intervals. To check accuracy the results were compared \v^ith 
another program which calculated the complex gamma function using 
a Pade-power approximation method and claimed 9-pl^ce accuracy. 
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The internal accuracy of the integration scheme was checked by 
comparing the 256 -point integrations to 512 -point integrations 
on the same intervals. This comparison indicated accuracy to 
about 13 places and these results were consistent with those found 
by the comparison program. In addition the numerical integrations 
seemed to take somewhat less time than was required by the power 
series method and was not subject to slow convergence as the ima- 



ginary part of the argument became large. 



COMPUTER PROGRAM 



PROGRAM PHASE MK IV SEPT 1972 

implicit REAM'S (A-H,0-Z) 

INTEGER RO 

EXTERNAL C GARG , F AC T , DC OS , D S I N , F QN , FC N 

STORAGE FOR SAMPLE POINTS AND WEIGHTS FOR 256 AND 512 POINT 
GAUSS-LEGEf.DRE INTEGRATIONS. 

C0MM0N/L256/X256(128) , A256 (128) 

COMMON/ L512/X5 12 ( 2 56 ) , A5 12 ( 256 ) 

INTEGRATION INTERVALS FOR GAMMA FUNCTION INTEGRATIONS 
COMMON/CD/DPC( 15) ,DMC( 15) 



MISCELLANEOUS CONSTANTS USED IN VARIOUS SUBROUTINES 
COMMON/ CONST/ PHI ,PH2, FRl , FR2, ASCL , P 1 , CA, CT , AA 



PASSES ALPHA AND AND ETA TO THE INTEGRATION AND MATRIX ROU- 
TINES 

COMMON /TRM/AL= A, ETA 



PASSES K TO THE INTEGRATION ROUTINE 
CUMMON/WAVP/WAVN 



PASSES NUCLEAR CHARGE 
CONMON/ATNO/Z 



DUMMY VARIABLES. THIS COMMON USED IN INTEGRATIONS FOR 
RHO = - 1 

COMMON/MCNE/Dl ,D2,D3 



USED IN DETERMINING THE LIMITS OF THE RANGE OF THE INTEGRAND 
OVER WHICH TO FIND THE ZEROES 

COMMON /I NT /XSM , FF, BO ,XLM 



DUMMY VARIABLES. THIS COMMON USED IN INTEGRATIONS FOR 
RHO > -1 

COMMON/GEZRQ/OA,D5,D6, 110, 120 
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I 



INDICES FO^ MATRIX ELEMENTS 

COMMON/ INDEX/N I , N J , L I , L J , M I , MJ , I S 

VALUE OF EXPONENT OF X IN BOUND-FREE INTEGRALS 
COMMON/RHO/RO 

STORAGE FOR BOUND-FREE INTEGRALS 

DIMENSION AS(l^f52) ,3S{ 1452) , AC (1452), BC( 1452) 

STORAGE FOR THE VALJES OF THE EXPONENTS OF Rl, R2 AND R12 IN 
THE HYLLERAAS FUNCTIONS. LOCATION OF THE REGIONS OF ALPHA 
TO BE SKIPPED IN THE CALCULATION 

DIMENSION N(95) ,L(95),M(95), ISTP(IO) 

DATA INPUTS FOR EXACT ENERGY CALCULATIONS 
DIMENSION X( 5) ,Y( 5) 

C DIMENSION STATEMENTS FOR PHASE SINGLET, N=7. 

STORAGE FOR BOUND-BOUND INTEGRALS 
DIMENSION A ( 5508) , B( 5508) 

STORAGE FOR MATRIX ELEMENTS AND EIGENVECTOR ELEMENTS 
DIME NS I ON H( 73 ,70) , S (7 J, 70 ) , C( 70 ,70 ) 

STORAGE FOR EIGENVALUES 

DIMENSION WAV( 70) , ENERG( 70) 

STORAGE FOR UPPER TRIANGULAR ELEMENTS OF H AND S AND FOR THE 
INVERSE OF THE TRIANGULAR FACTOR OF S 

DIMENSION TH (2485) , TS( 2485 ) 

EQUIVALENCE ( H ( 1 ) , AS ( 1 ) ) , (H (2451 ) , BS( 1 ) ) 

EQUIVALENCE ( S( 1) ,AC( 1) ) ,(S(2451) ,BC(D) 

C 

DATA FOR THE HYLLERAAS FUNCTIONS AND THE NUMERICAL INTEG- 
RATIONS ARE STORED IN A SEQUENTIAL DATA SET 

READ(8 , 1000 )N, L, M 

READ(8,1002) DPC,DMC 

READ! 8, 1 0 01 ) (X512( II), A512( II), 1 1 = 1,256) 

READ( 8, 1001) ( X 25 6( KK) , A256( KK) ,KK=1 , 12 8) 

REMAINDER OF DATA FOR PRESENT CALCULATION IS READ IN FROM 
PUN'CHED CARDS SUBMITTED WITH THE MAIN PROGRAM. THIS DATA 
REMAKES FIXED FOR THE ENTIRE PROGRAM: INDEX FOR NUMBER OF 

INTEGRATIONS, TOTAL SPIN, BASIS SET SIZE DETERMINED BY THIS 
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IfsDEX, NUMBER OF ENERGIES TO BE INVESTIGATED (SET THIS EQUAL 
TO ONE IF THE ENTIRE ELASTIC RANGE IS TO BE COVERED), 

NUCLEAR CHARGE 

READ! 5 ,1021) NSET, I S , HATRI X , I NO , Z 
IF( IS .EO.O ) GO TO 88 

IF TOTAL SPIN IS ONE (TRIPLET) SCREEN OUT THOSE FUNCTIONS 
NOT NEEDED FOR THE CALCULATION 

K = 0 

DO 100 KK=1,95 

IF(N(KK) .EQ.L( KK ) ) GO TO 100 
K = K + 1 
NIK) = N(KK) 

L(K) = L(KK) 

H(K) = M(KK) 

IF(K. EO. MATRIX ) KK = 95 
100 CONTINUE 

SET VARIOUS CONSTANTS WHICH DETERMINE NEEDED INTEGRALS 
88 MAX = 2"NslSET + 4 
I MAX = NSET 5 
MAXI = NSET + 3 
LIM = MATRIX^MATRI X 
I SIZE = MATRIX^ (MATRIX + l)/2 

FIND THE REQUIRED BOUND-BOUND INTEGRALS 
CALL TBL1(MAX,A,B) 

SET CONSTANTS NEEDED FOR THE NUMERICAL INTEGRATIONS 
PI = 0.314159265358979301 
XSM = 0.5D-2 
FF = 0.3D-1 
BO = 0.16D2 
XLM = 0.45502 

FIRST LOOP. READ DATA FOR THE RANGE OF ALPHA AND ENERGY TO 
BE INVESTIGATED. 

DO 130 IW=1,IN0 

READ (5 , 1008 ) I AL 0 , I AH I , I A I NT , NUM , FCT , PL 0 , P H I 
VALUES OF ALPHA TO BE SKIPPED 
REA0(5,1005) ISTP 
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USED TO LIMIT SEARCH FOR CLOSEST EIGENVALUES IN EXACT ENERGV 
CALCULATIONS. IF ND EIGENVALUE IS FOUl'JD WITHIN THIS EX- 
TENDED RANGE THE CALCULATION IS STOPPED FOR THAT ENERGY 

PHIEXT = PHI + O.lD-1 

PLOEXT = PLO - O.lD-1 

SECOND LOOP. COVERS DESIRED RANGE OF ALPHA FOR DESIGNATED 
ENERGY 

102 DO 131 I A=I ALO, I AHI , I AINT 

IF THE LOW VALUE OF K DESIRED IS ZERO, THE ENTIRE ELASTIC 
SCATTERING REGION WILL BE INVESTIGATED. SKIP THE EXACT 
ENERGY CALCULATION 

IFIPLD.EO. 3.3D3) GO TO 50 

READ THE VALUES OF ALPHA AND K FOR THE FIVE POINTS ON THE 
TRAJECTORY CLOSEST TO THE DESIRED VALUE OF K 

RE AD (5, 1316) (Y(I),X(I) ,1 = 1,5) 

XINT = PLO + 0.5D-6 

USING A CUBIC INTERPOLATING POLYNOMIAL, FIND AN ESTIMATE OF 
THE VALUE OF ALPHA REQUIRED TO GIVE THE DESIRED VALUE OF K 

70 CALL SPLIN1(X,Y,5,XINT,YINT) 

ALFA = YINT 

GO TO 51 

CHECK TO SEE IF A SEGMENT OF THE ALPHA RANGE IS TO BE SKIP- 
PED. NOT USED WHEN THE EXACT ENERGY CALCULATION IS USED 

50 DO 103 I=1,NUM,2 

IF ( I A. EO. I STP ( I ) ) IA=ISTP(1+1) 

103 CONTINUE 

ALFA = FCT*DFLOAT{ I A) 

EVALUATE THE H AND S MATRIX ELEMENTS 

51 DO 16 J = l, MATRIX 
J J = J=^ ( J - 1) /2 
DO 10 I =1 , J 

N I = N ( I ) 

NJ = N(J) 

LI = L( I ) 

LJ = L(J) 

M I = K ( I ) 

MJ = M(J) 
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CALL ELEM(HH,SS,A,B,MAX) 



I J 


= I 


+ JJ 


HI I 


f J ) 


= HH 


THI I J) 


= HH 


S I I 


, J ) 


= SS 


10 TSIIJ) 


= SS 



16 CONTINUE 

BEFORE FACTORING S, SCALE IT SO THE LARGEST ELEMENT IS < 10. 
F I S THE SCALE FACTOR 

CALL SCLITS, IS IZE, F) 

WRITE ( 6. 1006) F 

CALL SMPYI TS, F, ISI ZE) 

FACTOR S INTO ITS TRIANGULAR ELEMENTS AND INVERT THE UPPER 
TRIANGULAR FACTOR. RETURN IT IN TS 

CALL INVERT(TS,H,S,C,MATRIX, ISIZE) 

FORM U**(-T)*H*U**(-1). RETURN IT IN TH 
CALL TMPRDITH.TS.C, MATRIX, F) 

EXPAND TH TO FULL SIZE AND STORE IN H. PUT TS IN UPPER 
TRIA^;GLE OF S AND SAVE FOR FUTURE USE 

DO 2D J=l, MATRIX 

JJ = J'^IJ - l)/2 

DO 20 1=1, J 

I J = I + J J 

HT = H( I , J ) 

HI I , J) = THII J ) 

HI J , I) = THI IJ ) 

THI I J) = HT 

SI J, I) = O.ODD 

20 SII,J) = TSIIJ) 

TRACE = O.ODO 

DO 11 1=1, MATRIX 

11 TRACE = TRACE + HI I , I ) 

WRITE! 6 , 1013 ) TRACE 

FIND EIGENVALUES AND EIGENVECTORS OF TRANSFORMED H MATRIX 
CALL DJACVTIH,MATRIX,1 , E NE RG , C , MAT R I X ) 
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TRACE 



0.0 DO 



FIND SUM OF EIGENVALUES AND COMPARE WITH TRACE OF H PREV- 
IOUSLY FOUND 

DO 17 1=1, MATRIX 

17 TRACE = TRACE 4- ENERG(I) 

WRITEI6, 1015 ) TRACE 

SCALE EIGENVECTORS SO LARGEST ELEMENT IS < 10 
CALL SCL(C ,LIM,F) 

CALL SMPYIC, F, LIM) 

MULTIPLY THE ORTHOGONAL EIGENVECTORS BY U**(-l) TO FIND 
THE EIGENVECTORS OF THE ORIGINAL PROBLEM 

CALL MPRD(S,C, H, MATRIX) 

WRITE! 6,3002) ALFA 

WRITE (6,3004 ) ENERG 

IK = 0 

THIRD LOOP. EXAf'.INE THE EIGENVALUES TO FIND THOSE IN THE 
DESIRED ENERGY RANGE 

DO 15 I E=l, MATRIX 

WAV! IE) = O.ODO 

E = ENERG! IE ) 

ROOT =0.2D1=:E + Z 

IFIROOT .LE. O.ODO) ROOT = O.ODO 

PINT = DSQRT! ROOT) 

IF! (PINT.LE.PLO) .OR.IPINT.GT .PHI ) )G0 TO 21 

INDEX GIVING NUMBER OF EIGENVECTORS IN THE DESIRED ENERGY 
RANGE 

IK = IK + 1 
WAV! IK) = PINT 

EIGENVALUES IN DESIRED RANGE RELOCATED IN FIRST ADDRESSES OF 
ENIGENVALUE VECTOR 

ENERGIIK) = E 

EIGENVECTOR CORRESPONDING TO ABOVE EIGENVALUE RELOCATED IN 
APPROPRIATE COLUMN OF EIGENVECTOR MATRIX 

DO 12 K=l, MATRIX 

12 C! K, IK) = C !K, I E ) 



110 



IN EXACT EX’ERGY CALCULATION ASSUME FIRST VALUE WITHIN 
DESIRED LIMITS THE THE ONE LOOKED FOR 

IF(PLO.GT.O.ODO) IE = MATRIX 

GO TO 15 

IF EXACT ENERGY FEATURE IS NOT BEING USED SKIP THIS SECTION 
21 IF( PL3. EO. D.ODO) GO TO 15 

IF EIGENVALUE BEING EXAMINED IS OUTSIDE THE EXTENDED RANGE 
GO ON TO THE NEXT ONE 

I F ( ( PI NT. LT. PLOEXT ) .OR. {PINT . GT .PHIEXT ) )GO TO 15 

IF EIGENVALUE IS OUTSIDE THE DESIRED RANGE BUT WITHIN THE 
EXTENDED ONE, PRINT THE VALUE AND REPLACE THE OUTERMOST 
POINT ON THE SEGMENT OF THE TRAJECTORY UNDER EXAMINATION BY 
THIS VALUE. THEN GO BACK AND RECALCULATE AN IMPROVED 
ESTIMATE OF THE REQUIRED ALPHA 

WRITE(6.1014) IE, PINT 

KNUM = 0 

IF(PINT.LT.X(1 ) ) KNUM = 1 

IFCPINT .GT .X(5 ) ) KNUM = 5 

IF ( (KNUM.EO. 1 ) .OR. ( KNUM. EQ.5) ) GO TO 65 

DO 60 1=2,5 

IF((PINT.LT.X( I ) ).AND.(Pm.GT.X( I-IJ ) ) KNUM = I 

60 CONTINUE 

I F ( (KNUM.LT .1 ) .OR. (KNUM.GT .5 ) ) GO TO 15 
IF( { PINT-X(l) ) .GE. ( X( 5)-PINT) ) GO TO 61 
DO 62 1=1,4 
IL = 5 - I 

IF( IL.GE.KNUM) X(IL+1) = X(IU 

I F( IL .GE .KNUM) Y(IL + 1) = Y(IL) 

62 CONTINUE 
X(KNUM) = PINT 
Y(KNUM) = ALFA 
GO TO 70 

61 KNUM = KNUM - 1 
DO 63 1=2, 5 

IF(I .LE.KNUM) X(I-1 ) = X(I ) 

IF ( I .LE.KNUM) Y(I-1 ) = Y( I ) 

63 CONTINUE 

65 X(KNUM) = PINT 
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Y(KNUM) = ALFA' 

GO TO 70 
15 COrJTINUE 

IF NO EIGENVALUES WERE FOUND IN THE DESIRED RANGE, GO TO THE 
NEXT VALUE OF ALPHA 

I F ( I K. EO.O) GO TO 101 

WRITE(6,3005) ALFA 

WRITE! 6,30 04) (WAV( LP) ,LP=1 ,IK) 

FOURTH LOOP. BEGIN THE PHASE SHIFT CALCULATION 
DO 113 NN=1,IK 
E = ENERGI NN) 

WAVN = WAV(NN) 

FIRST FIND THE COULOMB PHASE SHIFT, ETA 
AA = (Z - O.IOD/WAVN 

IF Z = 1 SKIP THIS PART AND GO ON TO THE BOUND-FREE INTEG- 
RATIONS 

IFIAA. EO.O.ODO ) GO TO 301 

ASO = AA^AA 

CA =ASO - 0. 11D2 

CT = 0.6D1*(AS0 - O.IDD/AA 

CALL INT256(CSARG,GR) 

CA = CT 

CT= 0.1 1D2 -ASO 
CALL I NT256( C3 ARG ,GI ) 

GR = -GR 

ETA = DATAN2(G1 ,GR) 

GO TO 302 

301 ETA = O.ODO 

WRITE (6 ,1020) MATR IX, MATRIX, ALFA, E, WAVN 
GO TO 303 

302 WRITE! 6 , 1004) MATRI X , MAT RI X , ALF A , E , WAVN , ET A, Z 

FIND THE REQUIRED BOUND-FREE INTEGRALS 

303 CALL TBL2! IMAX,MAXI ,AS,AC,BS,BC ) 

STERM = O.ODO 

CTERM = O.ODO 
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FCRM THE BOUND-FREE MATRIX ELEMENTS 
DO 300 JM=1, MATRIX 
CC = CUM.NN) 

NJ = N( JM) 

LJ = L(JM) 

MJ = M( JH) 

CALL PHASE (AS, AC , BS , BC , MAX I , HS , HC , E ) 

HSC = HS=5^CC 
STERM = STERM + HSC 
HCC = HC'"*CC 
CTERM = CTERM + HCC 
300 CONTINUE 

EVALUATE THE PHASE SHIFT 

TDLTA = -STERM/CTERM 
TDLTAK = TDLTA/WAVN 
DLTA = DATAN(TDLTA) 

IF ( IS. EO.O) GO TO 91 
WRITE (6, 1011 ) 

GO TO 92 

91 WRITE(6, 1012 ) 

92 WRITE! o, 1007) ALFA , WAVN , TDLTA , TDLTAK ,D LTA 
WRITE (7 , 1100) ALFA, WAVN, TDLTA, TDLTAK, DLTA 
IFIDLTA.GT .O.ODO) GO TO 113 

IF THE PHASE SHIFT IS NEGATIVE, ARBITRARILY ADD PI TO IT AND 
PRINT AN ADDITIONAL DATA CARD WITH THE NEW PHASE SHIFT 

DLTAP = DLTA + PI 

WRI TE ( 7 , 11 01 ) ALFA, WAVN, DLTAP 

113 CONTINUE 

101 CONTINUE 

130 CONTINUE 

STOP 

1000 FORMAT (481 1 ,/, 47 11 ) 

1001 FORMAT ( 2D30. 16 ) 

1002 F0RMAT(8D9.3,/ ,7D9.3) 

1003 FORMAT (•-• ,5X, ' SOLUTIONS FOR ALPHA = • , F 1 1 . 6, / , ' - • ) 
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1004 FORMAT {' 1 lOX, ' ELECTRON - HELIUM + SCATTERING 
IPHASE SHIFT FOR', 13,' X ',12,' HAMILTONIAN MATRIX.',/, 
2 , 34X, 'ALPHA = ' , F 1 0. 7 , / , ' 0 ' ♦ 24X , • T OTA L ENERGY, E =', 

3024. 16, / , ' O' ,25X, ' WAVE NUMBER, K = ' , 024. 16 , / , ' O' , 1 5X , ' 
4C0ULCMB PHASE SHIFT, ETA = ' , 024 . 16 , / , ' 3 ' , 2 OX , NUCL E AR 
5CHARGE", Z =',F9.5) 

1005 FORNATdOIB) 

1006 FORMAT (' 1 ', lOX ,' F = ',010.3,/,'-') 

1007 FORMAT! ' O' ,5X, ' ALPHA = ' , F9 . 6 , 4 X , ' K = ' , 02 4 . 16 , 5 X , ' 
ITANGENT ( PHASE SHIFT) = ' , 024 . 16 , / , ' 0 ‘ , 65X , ' TAN( OE LTA ) /K 
2=' ,O24.16,/,'0' ,66X,'PHASE SHIFT =',024.16) 

1008 FORMAT ( 8X , 2 1 8 , 2 1 3 , F 1 0. 3 , 5X , 2F 1 5. 9) 

1011 FORMAT! ' O' ,/,'-' ,10X, ' TRI PLET PHASE SHIFTS.') 

1012 FORMAT ! 'O' ,/,'-' ,10X, 'SINGLET PHASE SHIFTS.') 

1013 FORMAT! lOX, ' TRACE OF ! U** ! -T ) ) U=!'=«' ! -1 ) ) =',024.16) 

1014 FORMAT! '-' ,10X, ' K! ', 12 ,' ) =',024.16,/,'-') 

1015 FORMAT ! 22X, ' SUM OF EIGENVALUES =',024.16,/) 

1016 F0RMAT!F11.7,F 12.7) 

1017 FORMAT! '-' ,10X, 15, ' CAROS PUNCHEO FOR THIS JOB.') 

1020 FORMAT! ' 1' ,/,'-' ,10X, ' ELECTRON - HYOROGEN SCATTERING 
IPHASE SHIFT FOR', 13,' X ',12,' HAMILTONIAN MATRIX.',/, 
2'-' ,34X, 'ALPHA = ' , F 9. 6 , / , ' 0 ' , 24 X , ' TOT AL ENERGY, E = ', 

3023. 16, /,' 0 ', 25X, ' WAVE NUMBER, K =',023.16) 

1021 FORMAT! 212 ,214 ,08. 3) 

1100 FORMAT ! Fll .7 ,F 12 .7 , 3X, 3018 . 10) 

1101 F0RMAT!Fli.7,F12.7,39X ,018. 10) 

3032 FORMAT !'-' ,5X, ' E IGENVALUES OF H FOR ALPHA = ',F12.7,/) 

3004 F0RM.AT!4X, 5025.16) 

3005 FORMAT! /,5X, ’ THE VALUES OF K FOR ALPHA =',F12.7,' 

1 ARE: ' , /, ' O' ) 

ENO 
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SUBROUTINE T B_ I { MA X , A , B I 



I 



TBLl GENERATES THE TABLE OF BOUND-BOUND INTEGRALS 
IMPLICIT REAL<=8 (A-H,0-Z) 

DIMENSION A( 1) ,B(1) 

THE INTEGRALS ARE STORED IN A ONE DIMENSIONAL ARRAY CORRES- 
PONDING TO THE VARIOUS VALUES OF NUt LAMDA AND MU 

MAX2 = MAX*MAX 

LIM = { MAX - 1 )*MAX2 

INITIALIZE THE A AND B ARRAYS 
DO 1 1=1, LIM 
A( I ) = . ODO 

1 B(I) = .ODD 
NMAX = MAX - 1 

EVALUATE THE INTEGRALS FOR MU = 2 
DO 2 NN=1,NMAX 
DO 2 LL=2,MAX 

IF{ (NN+LL.LT.5 ) .OR. (NN+LL.GT.MAX+2) J GO TO 2 
IN = NN + MAX^ILL - 1) + 2=?=MAX2 
CALL SUMINN.LL ,AA) 

A (IN) = AA 

2 CONTINUE 

USE THE RECURSION RELATIONS TO FIND THE INTEGRALS FOR MU = 
DO 3 NN=2,MAX 
DO 3 LL=2,MAX 

IF{ (NN+LL.GT.MAX+3 ) .OR. (NfI+LL.LT.6) ) GO TO 3 
IN = NN + MAX^ILL - 1) + MAX 2 
INP = IN + MAX 2 - 1 
A(IN) = A (INP) 

IF( ( LL. LT. 3) .OR. ( NN. LT.3) .OR. (NN+LL.LT .6 ) .OR. ( LL.GE. 
IMAX) .OR. (NN.GE .MAX) ) GO TO 3 

INP = I N + MAX* ( MAX + 1) - 2 

B( IN) = A( INP) / .3D1 

3 CONTINUE 

MMAX = MAX - 1 
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NLMAX = MAX - 2 



USE THE RECURSION RELATIONS TO FIND THE INTEGRALS FOR MU > 2 
DO 5 MM=A,MMAX 
DO 5 NN=2, NLMAX 
DO 5 LL=2, NLMAX 

IF( (NN+LL.LT .5 ) .OR . (NN+LL+MM.GT.MAX+5) ) GO TO 5 
IN = NN + MAX*(LL - 1) + (MM - 1)*MAX2 
INN = IN - 2^( MAX2 - 1 ) 

INL = IN - 2=f'^1AX*( MAX - 1) 

INNL = IN - MAX^M2*MAX - 1) + 1 
A(IN) = A(INN) + A(INL) - .2D1*B(INNL) 

IF( ( MH.GT.MAX-3) .OR. ( NN+ LL . LT . 6 ) . OR. ( LL . L E .2 ) . OR . ( NN « 
lLE.2).OR.(LL.3T.MAX-3) .OR. (NN.GT.MAX-3) ) GO TO 5 

B(IN) = (B(INN) + B(INL) - . 2 Dl>!^ A ( I NNL ) ) =«< D FL OAT ( MM-3 ) / 
IDFLOAT ( MM+1 ) 

5 CONTINUE 

RETURN 

END 



SUBROUTINE SUM ( NU , L MDA , CUT ) 

SUM EVALUATES A( MU , L AMDA , 2 ) FOR GIVEN NU, LAMDA 
IMPLICIT REALM'S (A-H,0-2) 

EXTERNAL FACT 
A = O.ODO 
NN = NU - 1 ■ 

DO 10 I=1,LMDA 
I M2 = LMDA - I 
INI = NN + IN2 

10 A = A + FACT( I N1 )/( FACT( IN2)*{ 0. 2D1*=«=( I N1 + 1))) 
LL = LMDA - 1 

OUT = FACT (LL) FACT (NN) - A) 

RETURN 

END 



ll6 
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DOUBLE PRECISION FUNCTION FACT(N) 



FACT FORMS N FACTORIAL 

IMPLICIT REAL<--8 (A-H,0-Z) 
FACT = .IDl 
IF(N.EO.O) GO TO 20 
DO 10 I=1,N 
FI = I 

10 FACT = FACT^FI 
20 RETURN 
END 
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SUBROUTINE S PL I N1 ( X , Y , M , X I NT , Y I NT ) 



SPLIN PROVIDES INTERPOLATED VALUE OF THE ORDINATE USING 
"CUBIC SPLINE FITTING" 

USAGE: 

FIRST CALL TO SUBROUTINE: 

CALL SPLINl ( X, Y,M, XI NTtYINT) 

SJBSEC'JENT CALLS USING THE SAME DATA 
CALL SPLINN(X,Y,N,X1NT,YINT) 

DESCRIPTION OF PARAMETERS 

X: MONDTONICALLY INCREASING ABSCISSA ARRAY 
Y: ON5-FOR-ONE CORRESPONDING ORDINATE ARRAY 
M: NUMBER OF X AND Y VALUES SUPPLIED < OR = 
BOD. (MODIFIED IN THE PRESENT CONTEXT TO 
NO MORE THAN FIVE) 

XINT: VALUE OF ABSCISSA FOR WHICH CORRESPOND- 
ING ORDINATE IS TO BE INTERPOLATED (OR 
EXTRAPOLATED) 

YINT: INTERPOLATED (OR EXTRAPOLATED) ORDINATE 

VALUE 

MATHEMATICAL METHOD: 

UPON FIRST ENTRY TO SPLIN, A CALL TO SPLICO IS 
MADE TO DETERMINE THE COEFFICIENTS TO BE USED IN 
PERFORMING THE INTERPOLATIONS. SEARCH FOR BRACKET- 
ING ABSCISSA VALUES IS ALWAYS MADE FROM THE REFER- 
ENCE LAST USED IN INTERPOLATING. 

REFERENCE: 

PENNINGTON, R.H., "INTRODUCTORY COMPUTER METHODS 
AND NUMERIGAL ANALYSIS", MACMILLAN, NEW YORK, 1965 

IMPLICIT REALMS ( A-H ) , RE AL (0-Z) 

DIMENSION X(M) ,Y (M) ,C( 4,6) 

CALL SPLICO( X, Y,M,C) 



ENTRY SPLINN(X, Y,M , XINT, YINT) 

3 IF (XINT-X( 1 ) ) 70,1,2 

70 K=1 

GO TO 7 

1 YINT=Y(1) 

RETURN 

2 IF (XINT-X( K+1) )6,4,5 

4 YINT=Y(K+1) 

RETURN 

5 K=K+1 

IF(M-K) 71,71,3 

71 K=M-1 
GO TO 7 

6 IF( XINT-X( K) ) 13 ,12 , 11 
12 YINT=Y(K) 

RETURN 
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13 K=K-1 
GO TO 6 

7 PRINT 101, XI NT 

101 FORMAT! 8H0XINT = E18.9.32H, OUT OF RANGE FOR 
1 INTERPCLAT ION) 

11 YINT=(X(K+l)-XINT)*(C(l,K)>i'(X( K+ 1 ) -X I NT ) « <" 2+C { 3 , K ) ) 
YINT=YINT+ (XINT-X(K) )=!' {C(2,K)*(XINT-X(K) )=^*2+C(4,K) ) 
RETURN 
END 



SUBROUTINE SPL I CO ( X ,Y , M, C ) 

IMPLICIT REAL ’^8 (A-H), REALMS (0-Z) 

DIMENSION X(M) ,Y(M) , C ( A , 6 ) , D ( 6 ) , P ( 6 ) , E { 6 ) , A { 6 , 3 ) , B ( 6 ) , 
1Z( 6) 

MM=M-1 

DO 2 K=1,MM 

D( K) =X( K+1 )-X( K) 

P(K)=D(K)/6. 

2 E(K)=( Y(K+1)-Y(K) ) /D(K) 

DO 3 K=2,MM 

3 B(K)=E(K)-E(K- 1) 

A( 1 ,2) =-l.-D( 1 ) /D( 2) 

A(1 ,3) = D(1 )/D( 2) 

A( 2, 3) =P( 2)-P{ 1)=!‘A{1,3) 
A(2,2)=2.=!'{P(l)+P(2))-Pm=!=A(l,2) 

A(2,3)=A(2,3)/A(2,2) 

B( 2)=B( 2)/A( 2, 2) 

DO 4 K=3,MM 

A(K,2) = 2.*(PK-1) + P(K))-P(K-1)vA(K-1,3) 
B(K)=3(K)-P(K-1)=«=3(K-1 ) 

A(K,3)=P(K)/A(K,2) 

4 B(K)=B(K)/A(K,2) 

0=D(H-2) /D(M-1 ) 

AIM, 1 ) = 1 . + 0+A( M-2, 3) 

AIM, 2) =-0-A(M, 1)*A(M-1 ,3) 

BIM) = B( M-2 )- A( M, 1 ) - BIM-1 ) 
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Z(M)=B( M) /A(M, 2) 

Mt\l=K-2 
DD 6 1=1, MN 
K = M-I 

6 Z(K)=B(K)-A(K, 3)=!=Z(K+1) 

Z(l) =-A(l,2)»!^Z (2 )-A(l,3)’^Z(3) 
DO 7 K=1,MM 
0=l./( 6.^D(K) ) 

C( 1 ,K} = Z(K )'^Q 

C( 2,K) = Z(K + 1)^=Q 

C(3,K)=Y(K)/D(K)-Z(K)*P(K) 

7 C(4,K]=Y(K+1)/D(K)-Z(K+1)*P(K) 
RETURN 

END 
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SUBROUTINE EL E ''U HH , SS t A , B , MAX ) 



ELEM GENERATES THE 'MATRIX ELEMENTS <CHI ( I ) | H | CHI ( J) > AND 
<CHI (I) ICHI ( J)> FOR H AND S, RESPECTIVELY 

IMPLICIT REAL1=8 (A-H.Q-Z) 

INTEGER PWR 

COMMON/ TRM/ ALFA, ETA 

COMMON/ATNO/Z 

COMMON/ INDEX/NI , NJ , LI , LJ ,M1 ,MJ , IS 

D I MENS ION A( 1) , B ( 1 ) , NU ( 52 ) , MUI 5 2 ) ,LMDA( 5 2) ,COEF( 52) 

NIJ = NI + NJ 

LIJ = LI + LJ 

NLIJ = NI + LJ 

LNIJ = LI + NJ 

MU = MI + MJ 

FNl = DFLOATINIJ + 1) 

FLl = DFLCATdJ + 1) 

FN = DFLOAT (NJ ) 

FL = DFLOAT (LJ ) 

FM = DFLOAT(MJ) 

MAX2 = MAX*MAX 

SET THE VALUES OF NJ , LAMDA, MU AND THE COEFFICIENTS FOR 
EACH OF THE 52 TERMS OF H(I,J) AND THE 4 TERMS OF S(I,J) 

NU( 1 ) = NI J 2 

NU( 3) = NU(1) 

NU ( 4 ) = NU ( 1 ) 

NU( 6) = NU( 1) 

NU( 8) = NU(1) 

NU( 9) = NU( 1) 

NU(37) = NU(1) 

LMDA(IO) = NU( 1) 

LMDA( 11) = NU( 1 ) 

LMDA(14) = NU( 1) 

LMDA(15) = NU( 1) 

LMDA( 16) = NU( 1 ) 

LMDA( 18 ) = NU( 1 ) 
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LMDA( 42) 


= NU( 1) 


NU( 2) = 


NIJ + 1 


LMDA( 12 ) 


= NU( 2) 


NU(39) = 


NU(2) 


LMDA(44 ) 


= NJ( 2) 


NU( 5) = 


NIJ 3 


NU(38) = 


NU(5) 


NU(40) = 


NU( 5) 


LMOA( 13) 


= NU( 5) 


LMDA(41 ) 


= NU( 5 ) 


LMDA( 43 ) 


= NU( 5) 


c 

II 


NIJ 


LMDA( 17 ) 


= NJ( 7) 


NU(IO) = 


LI J 1- 2 


NU(12) = 


NU( ID ) 


NU(13) = 


NU( 10) 


NU(15) = 


NU( ID ) 


NU(17) = 


NU( 10) 


MU(18) = 


NU( 10) 


NU(41) = 


NU ( ID ) 


LMDA( 1) 


= NU( 10) 


LMDA( 2) 


= NU( ID ) 


LMDA( 5) 


= NU( 10 ) 


U',DA( 6) 


= NU( 10) 


LMDA( 7) 


= NU ( 1 0 ) 


LMDA( 9) 


= NU( 10) 


LMDA08) 


= NU( ID) 


NU(ll) = 


LIJ + 1 


L.MDA( 3) 


= NU( 11 ) 


NU(43) = 


NU ( 11 ) 


LKDA(40 ) 


= NU( 11 ) 


NU(14) = 


LIJ + 3 


LMDAl 4) 


= NU( 14 ) 


NU(42) = 


NU( 14) 


NU(44) = 


NUd^t ) 


LMDA(37 ) 


= NU( 14 ) 




I « 




LMDA(39) 


= NU(14) 


NU(16) 


= 


L IJ 


LMDA( ; 


8) 


= NU( 16) 


n:u( 19) 


= 


NLIJ + 2 


NU( 21 ) 


= 


NU( 19 ) 


NU(22) 


- 


NU( 19 ) 


NU(24) 


= 


NU( 19 ) 


NU( 26) 


= 


NU( 19) 


NU(27 ) 


= 


NU( 19 ) 


NU(45) 


= 


NU( 19 ) 


LMDA(28) 


= NU( 19 ) 


LKDA(29 ) 


= NU ( 1 9 ) 


LMDA( 32) 


= NU( 19) 


LMDA(33 ) 


= NU(19) 


LM0A(34 ) 


= NU( 19 ) 


LMDA( 36) 


= NU( 19 ) 


LMDA( 50 ) 


= NU ( 1 9 ) 


NU( 20) 


= 


NLIJ + 1 


NU(47) 


= 


NU (23 ) 


LMDA(30 ) 


= NU(20) 


LMDA( 52 ) 


= NU(20) 


hU (23 ) 


= 


NLIJ + 3 


NU( 46) 




NU(23) 


NU(48) 


= 


NU (23 ) 


LMDA( 31 ) 


= NJ( 23 ) 


LM0A(49) 


= NU(23) 


LM0A(51 ) 


= NU(23) 


NU(25) 


= 


NLIJ 


LMDA(35) 


= NU(25) 


NU(28) 


= 


LNIJ + 2 


NU( 30) 


= 


NU( 23 ) 


NU(31 ) 


= 


NU(23 ) 


NU(23) 


= 


NU(28) 


NU(35) 


= 


NU(23 ) 


NU(36) 


= 


NU(28 ) 


NU( 49) 


= 


NU( 23 ) 



LMDA(19 ) 


= NU128) 


LMDA( 20) 


= NU(28) 


LMDA(23) 


= MU (28) 


LM0A( 24) 


= NU( 28) 


LMDA( 25) 


= MU(28) 


LMDA(27 ) 


= NU( 28 ) 


LMDA( 46 ) 


= NU( 28) 


NU(29) = 


LNIJ + 1 


NU ( 5 1 ) = 


NU( 29 ) 


LMDA(21) 


= MU( 29) 


LMDA(48 ) 


= NU(29) 


NU(32) = 


LNIJ + 3 


LMDA(22) 


= MU (32) 


NU(50) = 


NU( 32 ) 


NU(52) = 


NU( 32) 


LM0A(45 ) 


= NU(32) 


LMDA( 47) 


= NU(32) 


NU(34) = 


LNIJ 


LMDA( 26 ) 


= NU( 34 ) 


MU( 1) = 


MU 4- 2 


MU( 2) = 


MU ( 1 ) 


MU( 3) = 


MU( 1) 


MU ( 7 ) = 


HU ( 1 ) 


MU( 8) = 


MU ( 1 ) 


MU(IO) = 


MU( 1) 


MU ( 1 1 ) = 


MU ( 1 ) 


MJ(12) = 


U ( 1 ) 


MU (16) = 


MU ( 1 ) 


MU (17) = 


MU ( 1 ) 


MU(19) = 


MU( 1) 


MU(20) = 


HU ( 1 ) 


MU (21) = 


MU( 1) 


MU ( 2 5 ) = 


MU(1) 


MU (26) = 


MU ( 1 ) 


MU(28) = 


MU( 1) 


MU (29) = 


MU ( 1 ) 



i 




MU( 30) 
MU (34) 
MU (35) 
MU( 4) 
MU( 5) 
MU( 6) 
MU( 13) 
MU (14) 
MU( 15) 
MU (22) 
MU( 23) 
HU (24) 
HU ( 3 1 ) 
MU( 32) 
HU (33) 
MU (37) 
MU(38) 
MU (39) 
MU( 40) 
MU (41) 
MU (42) 
MU(43) 
MU (44 ) 
MU( 45) 
MU ( 4 6 ) 
MU ( 47 ) 
MU(48) 
MU (49 ) 
MU( 50) 
/•'.U (51) 
MU( 52) 
MU( 9) 
HU (18) 
MU( 27) 
MU (3 6) 
COEF( 



= MU(1)- 
= MU(1) 

= MU( 1) 

= MI J 
= MU(4) 

= MU(4) 

= MU (4) 

= MU (4) 

= MU(4) 

= HU (4) 

= MU(4) 

= MU(4) 

= MU (4) 

= HU(4) 

= MU (4) 

= HU (4) 

= HU(4) 

= MU(4) 

= MU(4) 

= MU (4) 

= MU (4) 

= MU(4) 

= MU (4) 

= MU(4) 

= MU (4) 

= MU (4) 

= MU(4) 

= MU(4) 

= MU(4) 

= MU(4) 

= MU (4) 

= MU f 1 
= MU (9) 

= MU(9) 

= MU(9) 

1) = -0. 25D0*ALFA*AL FA 
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COEF( 10) 


= 


COEF ( 1) 


CO£F(19 ) 


= 


COEF( 1 ) 


COEF( 28) 


= 


C0EF( 1) 


COEF( 2) 




0.5 DO^-ALFA’i'FNl 


COEF( 12) 


= 


COEF( 2) 


COEF(21) 


= 


COEF( 2) 


COEF(29 ) 


= 


COEF(2) 


COEF( 3) 


= 


0.5D0*ALFA=5^FL1 


COEF (11) 


= 


COEF(3 ) 


COEF( 20 ) 


= 


COEF( 3) 


COEFOO) 


= 


COEF ( 3) 


C0^F( 4) 


= 


0.5D0"ALFA*FM 


C0EF( 5) 


= 


C0EF(4) 


COEF(13) 


= 


C0EF(4) 


COEF( 14) 


= 


COEF( 4) 


COEF( 22) 


= 


C0EF(4) 


COEF(23 ) 


= 


C0EF(4 ) 


COEF( 31 ) 


= 


COEF( 4) 


COEF (32) 


= 


C0EF(4 ) 


CCEF(37) 




-COEF (4) 


COEF (38) 


= 


-COEF (4) 


CQEF(41 ) 


= 


-COEF (4 ) 


C0EF(42) 


= 


-COEF(4) 


COEF (45) 


= 


-COEF (4 ) 


C0EF(46) 


= 


-COEF (4) 


COEF (49) 




-COEF (4) 


COEF(50 ) 


= 


-C0EF(4 ) 


COEF( 6) 


= 


-FM*( FM + FN + 


COEF (15) 


= 


C0EF(6) 


C0EF(24) 


= 


C0EF(6) 


C0EF( 33) 


= 


COEF (6) 


C0EF( 7) 


= 


-0. 5D0=<=FN'’^FN1 


C0EF( 17) 




COEF( 7) 


COEF(26) 


= 


COEF(7) 


COEF( 34 ) 


= 


COEF(7) 


COEF( 8) 




-0. 5D0*^FL*FL1 



- z 



- z 



FLl) 
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COEF( 16 ) 


= 


COEF( 8 ) 


COEF( 25) 


= 


COEF( 8) 


C0EFI35) 


= 


C0EF(8) 


COEF( 9) 


= 


0. IDl 


COEF( 18) 


= 


COEF (9) 


C0EF(27 ) 


= 


C0EF(9) 


C0EF( 36) 


= 


COEF( 9) 


C0EFI39) 


= 


FM^'FN 


COEF( 44) 


= 


COEF( 39) 


C0EF(48) 


= 


COE F( 39 ) 


C0EF(51 ) 


= 


C0EF(39) 


COEF( 40) 


= 


FM^FL 


COEF (43) 


= 


COEF (40 ) 


C0EF(47) 




C0EF(40) 


COEF{ 52) 


= 


COEF(40 ) 


MAX2 = MAX* 


MAX 


HH 


0, ODO 


IF( IS . EO 


.0 ) 


GO TO 1 



FOP TRIPLET SCATTERING THE SIGNS OF THE LAST 26 TERMS IN 
H ( I , J) MUST BE CHANGED 

DO 1 1=19,36 

1 COEF( I ) = -COEF ( I ) 

DO 2 1=45,52 

2 COEF( I ) = -COEF( I ) 

THE VALUES OF NU , LAMOA AND MU FOUND ABOVE MUST BE CONVERTED 
TO A SINGLE INDEX CORRESPONDING TO THE STORAGE MODE OF THE 
INTEGRALS 

10 DO 6 K=l,36 

IN = NU(K) + 1 + MAX^LMDA(K) + MU(K)*MAX2 
PWR = NU(K) + LMDA(K) + MU(K) 

6 HH = HH *■ COEF(K)^A( IN) /ALFA*>^PWR 

DO 7 K=37,52 

IN = NU(K) + 1 + MAX*LMDA(K) + MU ( K ) =<N-1AX2 
PWR = NU(K) + LMDA(K) + MU(K) 

7 HH = HH + COEFIK )*B( IN)/ALFA**PWR 

AFTER FORMAT lO.N OF N(I,J) THE SAME PROCESS MUST BE REPEATED 
FOR S( I . J) 
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SS = O.ODO 
DO 9 K=l,10,9 

IN = NU(K) + 1 + MAX*LMDA{K) + MU(K)*MAX2 
PWR = NU{K) + LMDA(K) + MUCK) 

9 SS = SS + AC n )/ALFA**PWR 
IF CIS. EQ. 0) SO TO 11 
DO 13 K=19,28,9 

IN = NUCK) + 1 + MAX*LMDACK) + MUCK)*MAX2 
PWR = NUCK) + LMDACK) + MUCK) 

13 SS = SS - AC IM )/ALFA**PWR 
GO TO 12 

11 DO 8 K=19, 28,9 

IN = NUCK) +1 + MAX^LHDACK) + MUCK)*MAX2 
PWR = NUCK) + LMDACK) + MUCK) 

8 SS = SS + AC IN) /ALFA*=i^PWR 

12 RETURN 
END 
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SUBROUTINE SCUS,N,F) 



SCL FINDS THE ORDER OF MAGNITUDE OF THE LARGEST (ABSOLUTE 
VALUE) ELEMENT IN A GIVEN MATRIX 

IMPLICIT REAL<^8 (A-H,0-Z) 

DIMENSION S(l) 

TOP = DABSISd ) ) 

DO 100 I=1,N 

T = DABS (S ( I ) ) 

IF(T.GT.TOP) TOP = T 

100 CONTINUE 

FF, = DLOG10(TOP ) 

IPWR = IDINT(FF) 

F = 0.1D2=f='=(-IPWR) 

RETURN 

END 



SUBROUTINE SM=Y(A,C,N) 

SMPY MULTIPLIES A MATRIX BY A SCALAR, RETURNING THE SCALED 
MATRIX IN THE SAME LOCATION AS THE ORIGINAL MATRIX 

IMPLICIT REALMS (A-H,0-Z) 

DIMENSION A(l) 

DO 1 I=1,N 

1 A( I ) = A{ I )=l=C 

RETURN 

END 
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SUBROUTINE IN\/ERT( TS,H , S ,C ,N , I SI ZE) 



INVERT TAKES THE SYMMETRIC, POSITIVE DEFINITE MATRIX S 
(STORED IN THE COMPRESSED MODE IN TS) AND, USING DSINV, FAC- 
TORS IT ir.TO ITS UPPER TRIANGULAR FACTOR AND INVERTS THIS 
FACTOR. THEN USING THIS RESULT AS A FIRST APPROXIMATION THE 
INVERSE- IS IMPROVED BY AN ITERATION METHOD 

IMPLICIT REALM'S (A-H,0-Z) 

REAL=*=4 EPS 

DIMENSION TS(1 ) ,H( 1 ) ,S ( 1) ,C{ 1) 

KIT = 0 
EPS = O.lE-7 

CALL DSINV(TS,C,N,EPS,IER) 

I F{ lER . EO.-l ) GO TO 16 
ISIG = 1 

U**(-l) (IN TS) AND U (IN C) ARE MULTIPLIED TO FORM A TRIAL 
IDENTITY ^HICH IS CHECKED FOR OFF-DIAGONAL ELEMENTS WHICH 
MAY BE > 10=’'* (-20), AND STORED IN THE LOWER TRIANGLE OF H 

DO 3 K=2,N 

KK = K - 1 

KM = K*(K - l)/2 

DO 2 1=1, KK 

IK = K + N* ( I - 1 ) 

CC = O.ODO 

DO 1 J=I,K 

I J = ( J - 1 ) /2 + I 

JK = KM + J 

1 CC = CC + C( I J )=*TS ( JK) 

IF(DABS(CC ) .GT .O.lD-20) ISIG = 0 

2 H( IK) = -CC 

3 CONTINUE 



IF THE INVERSE IS ALREADY GOOD ENOUGH NO IMPROVEMENT IS 
NEEDED AND RETURN TO THE MAIN PROGRAM CAN BE MADE 

IF( ISIG.EO.l) GO TO 15 

14 KIT = KIT + 1 

THE CORRECTION TERM IS FORMED BY MULTIPLYING U*=!=(-l) AND THE 
TRIAL INVERSE 

DC 6 K=2,N 
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KK = K - 1 
KM = (K - 1) / 2 
DO 5 1=1, KK 
IN = N=S=(I - 1) 

IK = K + IN 
11=1+1 
CC = 0.000 
DO 4 J= I I , K 
IJ = J + I N 
JK = KM + J 

4 CC = CC + H( I J )'^TS ( JK) 

5 S( IK) = CC 

6 CONTINUE 

THE IMPROVED INVERSE IS FORMED BY ADDING THE CORRECTION TERM 
TO U**(-l) 

DO 8 K=2,N 

KK = K - 1 

KM = K=;^ ( K - 1) /2 

DO 7 I =1 ,KK 

IK = KM + I 

KI = K + N+(l - 1) 

7 TS( I<) = TS( K ) + S(KI ) 

8 CONTINUE 
IS IG = 1 

THE NEW TRIAL INVERSE IS FORMED AND TESTED FOR ANY ELEMENTS 
> 10^«(-20). THIS VALUE IS PUT IN THE LOWER TRIANGLE OF S 

DO 11 K=2.N 

KK = K - 1 

DO 10 1=1, KK 

IN = N«( I - 1) 

IK = K + IN 

CC = O.ODO 

II = I +1 

I F( I I .GT .KK ) GO TO 10 
DO 9 J=I I , KK 
IJ = J + IN 
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JK = K + N*(J - 1) 

9 CC = CC + H{ IJ )*H{ JK) 

IF(D<\3S(CC) .GT .O.lD-20) ISIG = 0 

10 S( IK) = CC 

11 CONTINUE 

IF{ ISIG. EC. 1) GO TO 15 



IF THE INVERSE IS STILL NOT GOOD ENOUGH THE NEW TRIAL IN- 
VERSE IS PUT IN H AVID THE PROCESS IS REPEATED 

DO 13 K=2,N 

KK = K - 1 

DO 12 I=1,KK 

IK'= K + N«n - 1) 

12 H( IK) = S( IK) 

13 CONTINUE 
GO TO 14 

15 WRITE! 6, 1001) KIT 
RETURN 

16 WRITE (6, 1000) 

STOP 

1000 FORMAT ( lOX Sr OP . I ER =-l, AND S CANNOT BE FACTORED.’) 

1001 FORMAT! lOX, ' INVERSE FOUND AF TE R ' , I 2 , ' I TE RAT 1 ONS . ' ) 
END 



SUBROUT INE DSINV(A,C,N,EPS,IER) 



DSINV INVERTS A GIVEN SYMMETRIC POSITIVE DEFINITE MATRIX 
(MODIFIED IN THE PRESENT CONTEXT TO PROVIDE ONLY THE INVERSE 
OF THE UPPER TRIANGULAR FACTOR OF A) 

DESCRIPTION OF PARAMETERS: 

A: DOUBLE PRECISION UPPER TRIANGULAR PART OF 
GIVEN SYMMETRIC POSITIVE DEFINITE N BY N 
COEFFICIENT MATRIX. ON RETURN A CONTAINS 
THE RESULTANT UPPER TRIANGULAR MATRIX IN 
D0U3LE PRECISION 

C: MATRIX IN WHICH THE UPPER TRIANGULAR FACTOR 
OF A IS STORED FOR LATER USE IN THE CALLING 
PROGRAM 

N; the NUMBER OF ROWS (COLUMNS) IN THE GIVEN 
MATRIX 

EPS: SINGLE PRECISION INPUT CONSTANT WHICH IS 
USED AS A RELATIVE TOLERANCE FOR TEST ON 
LOSS OF SIGNIFICANCE 

lER: RESULTING ERROR PARAMETER CODED AS FOLLOWS 
IER=0 NO ERROR 

I£R=-1 NO RESULT BECAUSE OF WRONG INPUT 
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PARAMETER N OR BECAUSE SOME RAD- 
ICAMD IS NON-POSITIVE (MATRIX A 
IS NOT POSITIVE DEFINITE, POS- 
SIBLY DUE TO LOSS OF SIGNIFI- 
CANCE) 

IER=K WARNING WHICH INDICATES LOSS OF 
SIGNIFICANCE. THE RAOICAND 
FORMED AT FACTORIZATION STEP K+1 
WAS STILL POSITIVE BUT NO LONGER 
> |EPS=^A(K+1,K+1) 1 

DOUBLE PRECISION A, DIN, WORK, C 

DIMENSION A(l) ,C(1 ) 

FIND THE UPPER TRIANGULAR FACTOR OF A 
CALL DMFSD(A,N ,EPS,IER) 

IF(IER) 9,1,1 

INVERT UPPER TRIANGJLAR MATRIX T 

PREPARE INVERSION LOOP 

1 IPIV = N=i (N + 1 )/Z 
IND=IPIV 

SAVE T BY DUPLICATING IT IN C 
DO 7 1 = 1 , IPIV 

7 C( I ) = A( I ) 

ITITIALIZE INVERSION-LOOP 
DO 6 I=1.N 
DIN=l.DO/A ( IPI V) 

A( IPIV) =DI N 
MIN=N 
KEND=I -1 
LANF=N-KEND 
IF(KEND) 5,5,1 

2 J=IND 

DO A K=1,KEND 

INITIALIZE ROW-LOOP 
W0RK=0 .DO 
MIN=.MIN-1 
LHOR=I PIV 
LVER=J 

START INNER LOOP 
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DO 3 L=LANF,MIN 

LVER=LVER+1 

LHOR=LHOR+L 

3 WORK=WDRK-s-A(L(/ER)’!'A(LHOR) 
A( J) =-WORK«DIIM 

4 J=J-MIN 

5 I PI V=I P I V-MI N 

6 IND=IND-1 
9 RETURN 

END 



SUBROUTINE DMF SD ( A , N , E PS , I ER ) 

DMFSD FACTORS A GIVEN SYMMETRIC, POSITIVE DEFINITE MATRIX 
INTO ITS UPPER AND LOWER TRIANGULAR FACTORS (WHICH ARE THE 
TRANSPOSE OF EACH OTHER). THE SOLUTION IS DONE USING THE 
SQUARE ROOT METHOD 0^ CHOLESKY. 

DIMENSION A(l) 

DOUBLE PRECISION DPIV,DSUM,A 

TEST ON WRONG INPUT PARAMETER N 
IF(N-l) 12,1,1 
1 IER=0 

INITIALIZE DIAGONAL LOOP 
KPI V=0 
DO 11 K=1,N 
KPI V=KPI V+K 
IND=KP IV 
LEND=K-1 

CALCULATE TOLERANCE 

TOL=ABS ( EPS*SNGL ( A (KP IV) ) ) 

START FACTORIZATION LOOP OVER K-TH ROW 
DO 11 I=K,N 
DSUM=O.DO 
IF(LEND) 2,4,2 
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START INNER LOOP 

2 DO 3 L=1,LEND 
LANF=KPIV-L 
LIND=I ND-L 

3 OSUM=DSUM+A(LANF )*A(LIND) 

TRANSFORM ELEMENT A(IND) 

A DSUM=A( I ND)-DSUM 
IF(I-K) 10,5,10 

TEST FOR NEGATIVE PIVOT ELEMENT AND FOR LOSS OF SIGNIFI- 
CANCE 

5 IF(SNGL(DSUM)-TOL) 6,6,9 

6 If’iDSUM) 12,12,7 

7 IF(IER) 8,8,9 

8 IER=K-1 

COMPUTE PIVOT ELEMENT 

9 DPIV=DSORT(DSJM) 

A( KPI V) =DPI V 

DP IV=1 .DO/ DPI/ 

GO TO 11 

CALCULATE TERMS IN ROW 

10 A( IND) = DSUM*DP IV 

11 IND=1ND+I 
RETURN 

12 IER=-1 
RETURN 
END 
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SUBROUTINE T MP RO ( H , U , C , N , F ) 



TMPRD FCRNS THE MATRIX PRODUCT U** ( -T ) ( - 1 ) AND THEN 

MULTIPLIES IT BY THE SAME FACTOR AS WAS USED TO SCALE S 

IMPLICIT REAL^^S (A-H,D-Z) 

DIMENSION H(l) ,U(1) ,C(1) 

LIM = N* (N + 1 )/2 

DO 3 L=1,N 

DO 3 K=L,N 

LK = K’MK - l)/2 + L 
UHU = O.ODO 
DO 2 I =1 .L 
HU = 0.000 

I L = LY{ L - 1) /2 + I 
DO 1 J=1,K 

IJ = J*( J - 1) /2 + I 

IF( I .GT . J) IJ = I* ( I - 1 )/2 + J 

JK = K'l'IK - l)/2 J 

1 HU = HU + H( IJ)=^U( JK) 

2 UHJ = UHU + U( IL )-HU 

3 C(LK) = UHU 
DO A I =1 ,LIM 

4 H( I ) = C( I )=^F 
RETURN 

END 
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SUBROUTINE DJ«^CVT ( A ,N , NQYE S , E I VU , E I VR , NDI M ) 



DJACVT CALCULATES I DOUBLE PRECISION ALL THE EIGENVALUES 
OF A REAL, SYMMETRIC MATRIX BY USE OF THE JACOBI VARIABLE 
THRESHOLD METHOD. CORRESPONDING ORTHOGONAL EIGENVECTORS ARE 
ALSO CALCULATED, IF DESIRED. ORIGINALLY PROGRAMMED BY W. 
POOLE, COMPUTER CENTER, U. C., BERKELEY. 

DESCRIPTION OF PARAMETERS 

A: THE REAL-8 SYMMETRIC (SQUARE) MATRIX WHOSE 
EIGENVALUES AND EIGENVECTORS ARE SOUGHT. ROW 
DIMENSION MUST BE NDIM. 

N: ORDER OF A, 0<N< /=ND I M<1 61 
NOYES: INDICATOR SUPPLIED BY USER. 

N0YES=0 EIGENVECTORS ARE NOT WANTED 
NOYES=l EIGENVECTORS ARE WANTED 
EIVU: THE EIGENVALUES ARE RETURNED IN THIS VEC- 
TOR, WHICH MUST BE DIMENSIONED TO AT 
LEAST N IN CALLING PROGRAM. THEY ARE IN 
NO PARTICULAR ORDER. 

EIVR: THE EIGENVECTORS (IF REQUESTED) ARE RE- 
TURNED AS THE COLUMNS OF THIS MATRIX, THE 
ROW DIMENSION OF WHICH MUST BE NDIM IN 
THE CALLING PROGRAM. THE COLUMN DIMEN- 
SION MUST BE AT LEAST N IN THE CALLING 
PROGRAM. THE VECTOR EIVR(I,J), I =1 , N 
CORRESPONDS TO THE EIGENVALUE EIVU(J). 

IF EIGENVECTORS ARE NOT REQUESTED, AN 
UNDIMENSIONED VARIABLE SHOULD BE PUT IN 
THE POSITION OF EIVR IN THE CALL STATE- 
MENT. 

NDIM: THE ROW DIMENSION (FIRST DIMENSION) OF A 
AND EIVR IN THE DIMENSION STATEMENT OF 
THE CALLING PROGRAM. NDIM<161 

REMARKS 

1) MATRIX A IS DESTROYED DURING CALCULATION. 

2) ALL EIGENVALUES ARE CALCULATED EACH TIME DJACVT 
IS CALLED 

3) VARIOUS DIAGNOSTIC MESSAGES MAY BE SENT TO STAN- 
DARD OUTPUT MEDIUM BY DJACVT. THESE INDICATE 
IMPROPER AFGUMENTS AND ARE GENERALLY SELF- 
EXPLANATORY. 

4) LIMIT 0= 160 FOR N (AND NDMI ) COULD BE RAISED BY 

CHANGING STATEMENT 21 

REFERENCES 

VDN NEUMAN, J., AMD GOLDSTEIN, "THE JACOBI METHOD 

FOR REAL SYMMETRIC MATRICES", JOURNAL OF ACM,1, 
1959 

WILKINSON, J. H., "THE ALGEBRAIC EIGENVALUE PROBLEM 
OXFORD UNIVERSITY PRESS, 1965, P. 277 FF. 

IMPLICIT REALM'S! A-H.O-Z) 

DIMENSION A ( ND I M , ND I M ) , E I VU ( N DI M) , E I VR ( N D I M , N D I M ) 

IF(N-1 )20,23,2 1 

20 PRINT 22,N 

22 FORMAT!* N= ',13,' IS TOO SHALL. LIMIT IS 1. RETURN 
ITO CALLING ROUTINE FROM DJACVT.') 

RETURN 

23 PRINT 24, A ( 1, 1 ) 

24 FORMAT (24H IN DJACVT , MATRIX A = ,E14.6) 

RETURN 
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21 IF(N-160) 1 , 1 ,3 
3 PRINT 2,N 

2 FORMAT! • N= ',13,' IS TOO LARGE. LIMIT IS 160. RETURN 
ITO CALLING ROUTINE FROM DJACVT.'I 

RETURN 

1 IF (NOYES) 99, 102 ,99 
99 CONTINUE 

DO 101 J=1,N 
DO 100 1=1 ,N 

100 EIVR( I , J ) = 0.0 

101 EIVR{ J, J)=1.0 

102 AT0P=0. 

D0'll2 J=1,N 
DO 111 I=1,J 

IF(A( I, J )-A(J, I ) )90,103,90 
90 PRINT 106, N,N 

106 FORMATIIAH IN DJACVT ( A ( , I 3 , 1 H , , 1 3 , 3H ) ) , ) 

PRINT 108,1 , J, J, I 

108 FORMAT!' A ! ' , I 3 , ' , ' , 13 , ' ) AND A ! ' , I 3 , ' , ' , I 3 , ' ) WERE 
lUNEQUAL. THEY WERE REPLACED BY THEIR MEAN IN DJACVT') 

A! I , J) = .5=!=!A! I , J )+A!J, I ) ) 

A! J, I)=A! I ,J) 

103 CONTINUE 

IF! ATOP-DABS! A ! I , J ) ) ) 106, 1 11 , 1 1 1 

104 ATOP=DABS! A! I , J) ) 

111 CONTINUE 

112 EIVU! J )=A! J ,J) 

I F ! ATOP )1 09, 109, 113 

109 PRINT 110 

110 FORMAT !26H IN DJACVT, MATRIX A = 0 ) 

RETURN 

113 AVGF = DFL0AT!N^!N-1) )*. 55 
D=0. 0 

DO 114 JJ=2,N 
DO 114 11=2, JJ 
S = A! I I-l , J J )/ATOP 

114 D=S*S+D 
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DSTOP=( 1.D-Q6) 

THRSH =OSQRT ( D/ AVGF) * ATOP 

115 IFLAG=0 

DO 130 JC0L=2,N 
JC0L1=JC0L-1 
DO 130 IR0W=1,JC0L1 
AIJ=A( IROW, JCOL) 

IF( DABS( AI J) -THRSH) 130,130,117 

117 AII=A( IROW, IROW ) 

AJJ=A( JCOL , JCDL) 

S=AJJ-AI I 

IF(DABS( AIJ )-l . D-09=^DABS( S) ) 13 0, 130, 118 

118 IFLAG=1 

IF( l.D-10* DADS (AIJ )-DABS(S))116,119,119 

119 S=. 7071067811865475 
C = S 

GO TO 120 

116 T=AIJ/S 

S=O.25/DS0RT (0 .25-!-T*T) 

C=DS0RT( 0. 5+S) 

S=2 .*T*S/C 

120 DO 121 1=1 , IROW 
T = A( I , I ROW) 

U=A( I , JCOL ) 

A( I , IROW)=C'-T-S*U 

121 A( I , JCOL ) = S*T + C=?U 
12= IROW + 2 

IF ( I2-JCCL) 127 , 127 , 123 
127 CONTINUE 

DO 122 1=12, j;OL 
T = A( I-l , JCOL) 

U = A( IROW, I-l) 

A( I-l , JCOL )=S*U+C*T 

122 A( IROW, I-l ) = C*J-S-T 

123 A( JCOL , JCOL) =S*AI J+C*AJJ 

A( IROW, IROW)=C*A( IROW, I ROW ) - S>t^ ( C*A I J- S* A J J ) 
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DO 124 J=JCOL, ''J 
T = A( I ROW , J ) 

U=A( JCDLiJ ) 

A( IROW, J)=C-T-S*U 

124 A( JCCL, J ) = S=^=H-C=<=U 
IF(NOYES) 131,126,131 

131 CONTINUE 

DO 125 1=1, N 
T = EIVR{ I , I ROW) 

EIVRd, IROW)=C<^T-E IVR( I,JCOL)*S 

125 EIVR( I , JCOL)=S=f^T+EI VR(I ,JCOL)>!=C 

126 CONTINUE 
S=AIJ/ATOP 
D=D-S=?'S 

IF( D-D STOP ) 1260, 129, 129 
1263 0=0. 

DO 128 JJ=2,N 
DO 128 11=2, JJ 
S=A( I I-l , JJ )/ATOP 

128 D=S*S+D 
DSTOP=( 1 .D-06) *D 

129 THRSH=DSORT{ D/AVGF)^ATGP 

130 CONTINUE 

IF( IFLAG)115,134,115 
134 T=A( 1,1) 

A( 1 , 1) =E IVUIl ) 

EIVU ( 1 )=T 
DO 132 J=2,N 
T=A( J, J ) 

A{ J, J)=EIVU(J) 

El VU( J) =T 
DO 132 1=2, J 

132 A(I-1, J)=A( J,I-1) 

133 RETURN 
END 
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SUBROUTINE MPR D ( A , B , C, N ) 



MPRD FORMS THE MATRIX PRODUCT A*B. THE RESULT IS RETURNED 
IN B. C IS USED FOR TEMPORARY STORAGE. 

IMPLICIT REALM'S (A-H,0-Z) 

DIMENSION A{ 1) ,B(1 ) ,C( 1) 

DO 2 1=1 ,N 

DO 2 K=1,N 

IK = N*(K-1) + I 

CC = O.ODO 

DO 1 J=1,N 

I J. = N* { J-1 ) + I 

JK = N*(K-1) 4- J 

1 CC = CC + A( IJ ) =’.'B( JK) 

2 C(IK) = CC 
LIH = N=i=N 

DO 3 I =1 , LIH 

3 B( I ) = C( I ) 

RETURN 

END 



SUBROUTINE INr256( FCT, F) 



INT256 IS A 256-POIvJT GAUSS-LEGENDRE QUADRATURE USING AN EX- 
TERNALLY DEFINED FUNCTION. THE SAMPLE POINTS, WEIGHT VALUES 
AND INTEGRATION SEGMENTS ARE READ INTO THE CALLING PROGRAM 
FROM AN EXTERNAL SOJRCE 

IMPLICIT REALM'S (A-H,0-Z) 

COMMON/L256/X( 128) ,A(128) 

COMMON/CD/DPC( 15) , DMC ( 15) 

F = 0. ODD 

DO 100 1=1,128 

ZZZ = O.ODO 

XX = X( I ) 

AA = A( I ) 

TO MINIMIZE COMPUTATION TIME THE INTEGRAND IS COMPUTED AT 
THE SAME RELATIVE POINT IN ALL 15 INTEGRATION INTERVALS 
FIRST AND THEN THE SUM OF THESE VALUES IS MULTIPLIED 6Y THE 
WEIGHT FACTOR FOR THAT POINT 

DO 10 J=l,15 

THE DPC AND DMC ARE, RESPECTIVELY, (0+0/2 AND ID-C)/2, 

WHERE D AND C ARE THE UPPER AND LOWER INTEGRATION LIMITS 

OP = DPC(J) 

DM = DMC(J) 

CHANGE THE VARIABLE TO CONFORM TO THE -1,+1 INTERVAL RE- 
QUIRED IN GAUSS-LEGENDRE QUADRATURE 

P = XX^DM 

ZP = DP + P 

ZM = DP - P 

10 ZZZ = ZZZ + DM*(FCT(ZP) + FCT(ZM)) 

100 F = F + ZZZ-AA 
RETURN 
END 



DOUBLE PRECISION FUNCTION CGARG(Z) 



CGARG EVALUATES THE INTEGRAND AT THE DESIGNATED POINT FOR 
THE GAMMA FUNCTION INTEGRATIONS 

IMPLICIT REALt^e (A-H,0-Z) 
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COMMOrj/CONST/PHl ,PH2 »FR 1 , FR2 , SCL , Ti-.OP I ,FCT2,FCTl ,A 
ARG = A*DLOG(Z) 

CGARG = ( Z^-3 )^DEXP (-Z )*( FCT2*DS IN(ARG ) + FCTl*DCOS(AR 

RETURN 

END 
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SUBROUTINE TBL2( I NA X , M AX I , AS , AC , BS , BC ) 



TBL2 TABULATES THE REQUIRED BOUND-FREE INTEGRALS 
IMPLICIT REAL<'8 (A-H,0-Z) 

C0MM0N/C0NST/PH1,PH2,FR1,FR2, ASCL, PI ,CA,CT,AAA 

COMMON/WAVP/WAV 

COMMON/TRM/ALFA, ETA 

COMMON/ ATNO/Z 

COMMON/RHO/N 

DIMENSION AC(1 ) , AS(1),BC(1),BS( 1) 

DIMENSION Al( 144) , A2{ 144) 

DI ME NS I ON SI (13 ) ,C2 (13 ) ,SS (11 ) , CC( 11 ) 

DIMENSION X0(500) 

ALFA2 = ALFA -H Z 

CALCULATE THE FREQUENCY AND PHASE FACTORS IN THE SINUSOIDAL 
TERMS 

FRl = WAV/ALFA2 

PHI = AAA*DLCG( . 2D1*FR1 ) + ETA 

FR2 = . 2D1'M-!AV/ALF A 

PH2 = AAA=^DLOG( .2D1*FR2) + ETA 

ASCL = 0.5D0*ALFA/ALFA2 

WRITE (6 ,3000) PHI, FRl , PH2 , FR2 , ASCL 

IF Z=1 THE NUMERICAL INTEGRATIONS NEED NOT BE DONE 
IF(AAA.EO.O.ODO) GO TO 104 

FIND THE ZEROES OF THE INTEGRAND FOR II 
CALL ZERO( 1 , L, XO, IQ) 

EVALUATE THE II INTEGRALS 
DO 100 1 = 1, I MAX 
N = I - 2 

CALL CS INT ( 1 ,L , 10, XO, ZZ, ASCL ) 

100 SKI) = ZZ 

FIND THE ZEROES OF THE INTEGRAND FOR IS 
CALL ZER0(2,L, XO, IQ) 
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EVALUATE THE IS INTEGRALS 
DO 101 1=1, MAX I 
N = I - 2 

CALL CSINT(2,L,IQ,X0,ZZ,0.1D1) 

101 SS(I) = ZZ 

FIND THE ZEROES OF THE INTEGRAND FOR 13 
CALL ZERO! 3,L, XO, IQ) 

EVALUATE THE 13 INTEGRALS 
DO 102 1=1 , IMAX 
N = I - 2 

CALL CSINT (3,L, IQ,XO, ZZ, ASCL ) 

102 C2( I ) = ZZ 

FIND THE ZEROES OF THE INTEGRAND FOR IC 
CALL ZER0(4,L, XO, 10) 

EVALUATE THE IC INTEGRALS 
DO 103 1=1, MAXI 
N = I ~ 2 

CALL CS INT ( A,L , 10, XO, ZZ,0. IDl) 

103 CC( I ) = ZZ 
GO TO 105 

FOR Z=1 EVALUATE THE REQUIRED INTEGRALS EXACTLY 

104 CALL HINT! SI, C2, SS,CC, IMAX, MAXI ) 

105 WRITE (6, 2001 ) 

WRITEI6, 2000 )(J,S1(J),C2(J),J=1, IMAX) 
WRITE! 6,2002) 

WRITE(6,2000)(J,SS (J), CC ( J ) , J= 1 , MAX I ) 

MAXI I = MAXI 1 

INITIALIZE THE AS, AC, BS AND BC ARRAYS 
DO 199 J=l,1452 
AS(J) = .ODO 
BS(J) = .ODO 
AC(J) = .ODO 
199 BC(J) = .ODO 
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INITIALIZE THE A1 AND A2 ARRAYS 
DO 200 J=l,144 
AKJ) = .ODO 

200 A2(J) = .ODO 
KNMAX = MAXI 
LLMAX = MAXI 
NNM = NNMAX + 1 
LLM = 2*LLMAX - 1 
AP = 0.5D0-AL=A 
ALFAP = 0. IDl/ AP 

ALFAA = 0.2D1/ (0.2D1=!'Z + ALFA) 

USING THE IS AND IC INTEGRALS FOUND ABOVE EVALUATE AS(NU, 
LAMDA.2) AND AC( NU , L AMDA , 2 ) 

DO 201 NN=1, NNMAX 

SSN = SS(NN) 

CCN = CC(NN) 

AP = AP*ALFAP 

F = O.lDl 

ALFAAP = O.lDl 

DO 201 LL=1, LLMAX 

ALFA4P = ALFA4P'=ALF A4 

F = F=:<DFLOAT{ LL) 

IF ( (NN+LL.LT ) .OR . (NN + LL .GT . IMAX ) ) GO TO 201 
INA = NN + NNM*(LL + LLM) 

TAA = F^ALFA4P*AP 
AS(INA) =TAA<=SSN 
AC(INA) =TAA<'CCN 

201 CONTINUE 
NIMAX = MAXI I 
LI MAX = MAXI I 

I MAX I = I MAX 1 

USING THE II AND 13 INTEGRALS FOUND ABOVE EVALUATE AKNU, 
L AMDA, 2) AND A2 { NU , L AMDA , 2 ) 

DO 202 NN=1, NIMAX 

DO 202 LL=1 ,L1 MAX 

IF { (NN + LL.GT . IMAXI ) .OR . (NN+LL . LT. 5) ) GO TO 202 
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IN 12 = NN + N1MAX=!<(LL - 1) 

CALL A12SUM( NN, LL, AA, S 1 ) 

Al( IN12) = Ai\ 

CALL A12SUM(NN, LL, AA.C2) 

A2(IN12) = AA 

202 CONTINUE 
NNMAX = MAX I I 
LLMl = LLMAX - 1 

NLM = NNMAX^LLMAX - 1 

NNM = NNMAX - 1 

LNM = NNMAX* (-LM AX + 1) - 2 

i 

USING THE RECURSIOrj RELATIONS EVALUATE AS ( NU , L AM DA , 1 ) , 
AC(NU, LAMDA, 1 ) , BS ( 'J U , LAMDA , 1 ) AND BC ( NU , LAMDA , 1 ) 

DO 2 03 NN=2,NfJMAX 

DO 203 LL=1, LLMAX 

1F( (NNI + LL.GT.I MAXI ) .OR . (NN + LL. LT.5) ) CO TO 203 

INI = NN + NNMAX* (LL + LLMlI 

INAl = INI + NLM 

I NLA = LL + N1MAX*(NN - 1) 

INNA = LL + 1 + N1MAX*(NN - 2) 

AS(INl) = AS(INAl) + AKINLA) - AKINNA) 

ACIINl) = AC(INAl) + A2(INLA) - A2(1NNA) 

IF ( (NN .LT . 3) .OR . ( LL .LT .2) .OR . ( NN + LL. LT. 6) . OR. ( NN.GT. 
INNM) .OR. (LL.GT.LLM) ) GO TO 203 

INBl = INI + LNM 

INLB = LL - 1 + N1MAX*NN 

INNB = LL + 2 + N1MAX*-(NN - 3 ) 

BS(lNl) = ( AS(INBl) + Al(INLB) - A1 ( I NNB ) ) / . 3 D1 
BC(INl) = ( AC(INBl) + A2(INLB) - A2 ( I N NB ) ) / . 3D 1 

203 CONTINUE 



MMMAX 


= MAXI 


MAX 13 


= IMAX + 3 


MMM = 


MMMAX - 2 


NNl = 


2* (NNMAX*LLMAX 


NN2 = 


2*NNMAX*LLM1 


NM3 = 


NNMAX*LLM - 1 
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USING THE RECURSION RELATIONS EVALUATE AS ( NU , L AMDA, MU ) , 
AC(NU, LAMDAtMU), BS ( NU , L AMDA , MU ) AND BC ( NU , L AMDA , MU ) FOR 
MU > 2 

DO 204 MM=4,MMMAX 

FM31 = DFLOATIMM - 3I/DFL0ATIHM + 1) 

DO 204 NN=2,NNMAX 
DO 204 LL=1,LLMAX 

IF ( (NN+ LL+MM.GT .MAXI3 ) .OR. (NN+LL.LT.4 ) ) GO TO 234 
IN = NN + NNMAX*((LL - 1) + LLMAX*(MM - 1)) 

INN = IN - NNl 
INL = IN - NN2 
INNL = IN - N^J3 

AS(IN) = AS(INN) + AS(INL) - BS ( INNL ) * . 2D 1 
AC(IN) = AC(INN) + AC(INL) - BC ( INNL ) ^ . 2D 1 
IF( ( LL. LT. 2) .OR. (NN.LT.3) .OR. ( LL + MM.LT.6) .OR. (MM.GT. 





IHHM)) GO TO 


204 








BS( IN) = ( 


BS! I NN) 


+ BS!INL) - 


AS ! INNL )=?'.2D1 )=?'FM31 




BC(IN) = ( 


BC ! INN ) 


+ BC!INL) - 


AC! INNL )>•'. 2D1) 4FM31 


204 


CONT INUE 










RETURN 








2000 


FORMAT! 20X , 


12, 2D30. 


16) 




2331 


FORMAT! //// 


,3^X, ' 


SI' ,24X , ' 


C2 ' ,// ) 


2002 


FORMAT !//// 


,3^X, ' 


SS ' ,24X, ' 


CC ,//) 


3000 


FORMAT ! '-' , 


lOX , 'PHI 


= ' ,D23. 16 ,/ , 


' O' , lOX , ' FRl = ' , D23 



116,/, ' O' ,10X,' PH2 = ' , D23. 16,/ , • O' , lOX, ' FR2 = ',D23.16 
2,/,' 3' ,13X,'ASCL = ',022.16,/,'-') 

END 
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SUBROUTINE ZER 0( I C S , J , XO , I Q) 



ZERO FINDS ALL THE ZEROES OF SIN(A*X + B*LN(X) + C) OR 
COSIA^X + B*LN(X) + C) BETWEEN X=XSM AND X=SLM 

IMPLICIT REALM'S (A-H,0-Z) 

COMMON/CONST/PHl,PH2,FRl,FR2i ASCL.PI ,CA,CI ,B 
COMMON/INT/XSM.FF ,B0 ,XLM 
DIMENS ION XO(l ) 

X = XSM 

GO TO (21,22,21,22 ) , ICS 

21 A = FRl 
CC = PHI 
GO TO 1 

22 A = FR2 
CC = PH2 

1 ARG = A*X + BF'DLOG(X) + CC 

FIND THE SMALLEST (ALGEBRAICALLY) INTEGER MULTIPLE OF PI 
SUCH THAT Avx + B=^LM(X) + C > N’i'PI 

ARGP = ARG/PI 

ITHTA = IDINT(ARGP) - 1 

IT = IABS( ITHTA) 

DETERMINE WHETHER THIS ZERO STARTS A POSITIVE OR A NEGATIVE 
SWING OF THE FUNCTION AND RECORD THIS FACT 

10 = MODI IT, 2) 

IF( IQ.EO.O ) U = -1 

THTA = PI^i^DFLOATI I THTA) 

IF (ICS.GT.2) THTA = THTA - 0.5D0*PI 

FIND THE VALUE OF X CORRESPONDING TO THIS INITIAL VALUE OF 
THE ARGUMENT 

CALL XZERO(THFA,X, A,B,CC) 

X0(1) = X 

REPEAT THIS PROCESS UNTIL X > XLM OR UNTIL 500 ZEROES HAVE 
BEEN FOUND 

DO 2 J=2,500 

THTA = THTA + PI 

X=X + PI/(A+ B/X) 
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CALL XZEROCTHf A, X, A,B, CC) 
XO(J) = X 

IF (X.GT.XLM) G0T03 

2 CONTINUE 
J = 500 

3 CONTINUE 
RETURN 
END 



SUBROUTINE X ZE RO ( THTA , X , A , B » CC ) 

XZERO SOLVES THE EQUATION THETA = A*X + B>:aN(X) + C FOR X BY 
AN ITERATION PROCESS 

IMPLICIT REAL<'8 (A-H,0-Z) 

1 XI = X*(THTA - (CC + B«(DLOG(X) - 1)))/(X*A + B) 

IF THE LINEAR APPROXIMATION RESULTS IN A NEGATIVE ESTIMATE 
FOR X THE METHOD FAILS SO REPLACE THE ESTIMATE OF X BY ONE- 
TENTH OF THE INITIAL VALUE OF X AND REPEAT THE APPROXIMATION 
FINDING A NE.'. ESTIMATE. REPEAT THIS PROCESS UNTIL THE 
ESTIMATE BECOMES POSITIVE 

I F(X1 .LE.O .ODD ) XI = 0.1D0*X 

IF (DABSIX-Xl) .LT.O.lD-12) GO TO 2 

X = XI 

GO TO 1 

2 X = XI 
RETURN 
END 
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SUBROUTINE CSINT(ICS,M,IQ,XO,F,PI 



CSINT CARRIES OUT THE NUMERICAL INTEGRATION OF THE BOUND- 
FREE INTEGRALS FOR I > 1 

IMPLICIT REALMS (A-H.O-Z) 

EXTERNAL FQN, FCN , DS I N, DCOS 

COMMON/MONE/C, D,FFA 

COMMON/ I NT/XSM,FF, B0,XLM 

COMMON/GEZRO/FA, B, AMP , ICC, I QI 

COMMON/RHO/N 

DIMENSION STS( A ) ,T I ( 3) , AC( 4) , ADI 4) , AES( 4) , XO(M) 

ICC = ICS 
TS = O.ODO 
SINT = O.ODO 
AESI 1) = O.lDl 
DO 4 LL=1,3 

4 AESI LL+1 1 = AE SI LU + P 
FN = DFLOATIN) 

SINCE THE REGION WHERE THE INTEGRAND IS SIGNIFICANTLY DIFF- 
ERENT FROM ZERO IS DIFFERENT FOR DIFFERENT N, ADJUST THE 
POINT WHERE THE CAREFUL NUMERICAL INTEGRATION STARTS AND 
ENDS TO TAKE ADVAf.TAGE OF THIS 

XS = XSM 

IFIN.GT.l) XS = FF^IFN - D . 1 01 ) * I FN*=<^ 0 .6D D ) 

XL = 0.25D1'-FN + BO 
I FI N.LE. 1) GO TO 7 
DO 5 J=1 ,M 

IFIXOI J) .GE.XS) GO TO 6 

5 CONTINUE 

6 JJ = J - 1 

IFI JJ. EQ.O) JJ = 1 

lOI DETERMINES WHETHER THE FUNCTION IS STARTING A POSITIVE 
OR NEGATIVE SWING ON THE NEXT HALF-CYCLE 

IQI = M0DIJJ,2) 

IFIIQI.EO.l) IQI = IQ 

I FI IQI . EO. 0) I 01 = -IQ 

GO TO 8 
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7 JJ = 1 
lOI = 10 

8 D = X0( JJ ) 

IF (N.EO.-l) 30 TO 99 
DO 9 LL=1,4 

IFID.EO. 0.000) GO TO 100 
AD(LL) = DEXP(-AES (LL)*D)/D 
GO TO 9 

100 AD(LL) = O.ODO 

9 CONTINUE 

99 DO 10 J=1,M 
JI = M + 1 - J 
IF (X0( JI ) .LT. XL) GO TO 11 

10 CONTINUE 

11 IF( JJ .GE.J I ) J I = M - 1 
DO 17 K=JJ,JI 

C AND D REPRESENT SUCCESSIVE ZEROES OF THE TRUE INTEGRAND 
C = D 

D = X0(K+1) 

FA IS THE FREQUENCY OF THE SIN FUNCTION WITH ZEROES AT C AND 
D 

FA = 0. 3141592 653589793D1/ (D - C) 

FFA = FA 

DIO = DFLOATI I 01 ) 

IF THE TRUE INTEGRAND IS POSITIVE THE APPROXIMATE FUNCTION 
IS REDUCED IN AMPLITUDE BY 10? AND IF NEGATIVE, INCREASED 
BY 10? 

AMP = 0. IDl + 0. IDO - DIO 
IF(N.EO.-l) GD TO 14 
DO 12 LL=1,4 
AC(LL) = AD(LL) 

12 AD(LL) = DEXP( -AES ( LL ) *D )/D 

FIRST THE INTEGRAL 3F THE APPROXIMATE FUNCTION IS EVALUATED 
DO 13 LL=1,4 
AE = -A.ES( LL ) 

AEC = AC(LL) 
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AED = AD(LL) 

CALL CSSUM(D,C,FA,TERM,AE,AEC, AED) 

13 STS(LL) = TERM 

TS = TS + DIO>:'AMP=^-(STS( 1) - 0. 3D1*!' ( STS (2 ) 

1 - STS (3 ) ) - STS (4 ) ) 

GO TO 16 

14 UL = 0. IDl 

FOR N=-l THE APPROXIMATE FUNCTION MUST ALSO BE EVALUATED 
NUMERICALLY 

DO 15 LL = 1 ,3 

BL = UL 

UL = BL + P 

CALL D0G32 ( BL, UL ,FQN,TT) 

15 TI ( LL) = TT 

TS = TS - OIQ=?'AHP*FA>!=(TI ( 1) - 0.2D1«TI(2) + TK3)) 

NEXT EVALUATE THE INTEGRAL OF THE DIFFERENCE BETWEEN THE 
TRUE INTEGRAND AND THE APPROXIMATE ONE NUMERICALLY USING A 
32-POINT GAUSS-LEGENDRE QUADRATURE 

16 DD = DFLOATd + (1 - IOD/2) 

CC = DD -I 0.13 1 

B IS THE PHASE OF THE APPROXIMATE FUNCTION, SO CHOSEN TO 
MAKE THE FUNCTION BEHAVE PROPERLY 

B = FA*(DD'"D - CC>!=C) 

CALL D0G32 ( C,D, FCN, FG) 

SINT = SINT + FG 

17 lOI = -lOI 

FINALLY THE REGION -ROM 3 TO XS AND FROM XL TO 100 IS INT- 
FGPAT5D WITHOUT USING THE APPROXIMATE FUNCTION, UTILIZING A 
512-PCINT GAUSS-LEGENORE QUADRATURE 

IF ( ICS.GT.2) GO TO 18 

D = XO(JJ) 

IF (D.EO.O.ODO) GO TO 24 
CALL DC512( O.ODO,D,DSIN,ZI ) 

GO TO 25 

24 ZI = O.ODO 

25 C = X0(JI+1) 

IFIC.GE .0.6D2) GO TO 20 
CALL DQ512( C ,0. 6D2 , DSI N,UI ) 
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CL = 0.602 
GO TO 21 

20 CL = C 

UI = 0.000 

IF(C.GE. 0.103) GO TO 26 

21 CALL 0Q512(CL,0. 103, OSIN, UJ) 

GO TO 19 

18 0 = XO ( JJ ) 

IF(O.EQ. 0.000) GO TO 28 
CALL 00512(0. 000, 0,0C0S,Z1) 

GO TO 29 

28 Z1 = 0.000 

29 C = X0( J I+l ) 

IF(C.GE.O. 602) GO TO 22 
CALL 0Q512 (C,0 .602 , OCOS,UI ) 

CL = 0.602 

GO TO 23 

22 CL = C 

UI = 0.000 

IF(C.GE. 0.103) GO TO 26 

23 CALL OQ512(CL,0. 103,0C0S,UJ) 

GO TO 19 

26 UJ = 0.000 

THE FINAL INTEGRAL IS THE SUM OF ALL THE INOIVIOUAL INTEG- 
RATIONS 

19 F = SINT + TS + ZI + UI + UJ 
RETURN 

ENO 
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SUBROUTINE C S5 UM (D , C , A ,T ERM, AE , AEC , AED ) 



CSSUM FINOS THE INTEGRAL OF ( X'^=!=N ) ^ ( EX P { - AL PHA *X ) ) -S IN ( AX + B ) 
BETWEEN SUSCESSIVE ZEROES OF THE INTEGRAND 

IMPLICIT REAL<=8 {A-H,0-Z) 

COMMON /RHO/N 

LIM = N + 1 

DEN = -DSORT(A*A + AE*AE) 

PH = DATAN2I A, AE) 

S = O.ODO 
F = 0. IDl 



DI 


r= 


0. IDl 


DR 


= 


0 .IDl/DEN 


PHP 




PH^DFLC!^ T( N+2) 


FN 


= 


0. IDl/DFLOAT (LIM) 


AC 




AEC 


AD 


= 


AED 


DO 


10 


1=1, LIM 


F I 


= 


DFLCAK I ) 


FN 


= 


F N* F I 


DI 


= 


DI*DEM 


DR 


= 


DR* DEN 



F = F’i^FI 

PHP = PHP - PH 

AC = AC-C 

IF ( (AC. EQ.O. ODO) .AND. ( I .EQ.l ) ) AC = O.lDl 
AD = AD^O 

10 S = S DSIN(PHP)*FI*DR>^{AD + AO/F 
TERM = FN*S/DI 
RETURN 
END 
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SUBROUTINE DOS 32 ( X L , XU , FCT , Y ) 



D0G32 IS A 32-PDIVT GAUSS-LEGEN3RE QUADRATURE WHICH INTEG- 
RATES AN EXTERNALLY DEFINED FUNCTION FCT BETWEEN THE LIMITS 
XL AND XU 

DESCRIPTION OF PARAMETERS 

XL: LOWER BOUND OF THE INTERVAL 
XU: UPPER BOUND OF THE INTERVAL 

FCT: NAME OF THE EXTERNALLY DEFINED FUNCTION 
TO BE INTEGRATED 
Y: RESULTING INTEGRAL VALUE 

DOUBLE PRECISION XL , XU , Y , A , B , C , FCT 

DEFINING A AND B ENABLES THE VARIABLE TO BE CHANGED TO CON- 
FORM WITH THE -1,+1 INTERVAL REQUIRED OF THE GAUSS-L EGENDRE 
QUADRATURE 

A=..5D0*( XU+XL) 

B=XU-XL 

C=. 4986 3 193092 47A 07 6DO*B 

Y = . 35093 05 0 047 350483D-2’i^( FCT (A + 0 + FCT (A-C ) ) 

C= .49 280 57 557 726 3417D0*B 

Y = Y+.81371973S 54 52 835D-24 ( FCT I A+ C ) +FCT ( A-C ) ) 
C=.48238112779375322D0*B 

Y=Y+.l 26960326 5463 103 00-1=4-- (FCT (A^C) + FC T( A-C) ) 

C = . 4674 530 3796 8 86984 DC 

Y=Y-h. 17 1369314565107 170-1=!^ ( FCT (A 40 ■^FCT( a- C) ) 

C = . 4481 6 057788 3 026 06DO==B 

Y=Y4-.21 41 794901 1 113340 0-1=!' ( FCT { A4-C ) 4-FCT ( A-C ) ) 
C=.42468380686628499D0*B 

Y=Y4-. 2 5499 0296 31 186 088 0-1 « ( FCT { A4-C ) 4-FCT ( A-C ) ) 

C= .3972418979C39712GD0=’'B 

Y=Y4-. 29 342 0467 3926 7774 D-1 4= ( FCT ( A4-C) 4-FCT ( A-C) ) 

C= .366 09 10 5937 01448 4D0*B 

Y=Y4-. 32911 11 15 881809230-1-"= (FCT (A4-C)4-FCT( A-C) ) 

C=. 33 15221334651 076OD0*B 

Y = Y4-. 36 1728970 54 42 4253D- 1* I FCT { A4-C ) 4-FCT { A-C ) ) 

C=. 293 857878620581 16D0«B 

Y= Y4-. 390969473 93 535 1530-1=:-' { FCT ( A4-C ) 4-FCT ( A-C ) ) 

C= .25 34499 5446 61 14 70D0=!'B 

Y=Y4-. 4 16559621 1347 3378 D-l*( FCT (A4-C) 4-FCT (A-C) ) 

C= .21067 5 63806 53I767D0'-:'B 
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Y=Y+.43826 J465J220L9D6D-1* ( F CT ( A+C ) + FCT ( A- C ) ) 

C= .155934301 1^1063 82D0>:'B 

Y = Y+. 45 5 869393 47881 9420- I’i' ( FCT( A+C 1 + FCT ( A-C) ) 

C = . 1196436811? 5 068 54D0»'B 

Y = Y + . 469 22 1995 404 022830- !=:=( FCT ( A + C) + FCT ( A-C) ) 

C=. 72235980791 3982 50-1* B 

Y =Y+.47 8 193 60 03963 74300-1* ( FCT ( A+C) +FCT{ A-C) ) 

C=. 241 53832 8438691580-1*8 

Y= B* ( Y+. 482700 4425 7363 90 00-1 *( FCT ( A+C ) +FCT ( A-C ) ) ) 

RETURN 

ENO 



OOUBLE PRECISION FUNCTION FQN(Z) 

FON EVALUATES THE APPROXIMATE FUNCTION AT GIVEN Z FOR RH0=-1 
IMPLICIT REAL*8 (A-H.O-Z) 

COMMON/ MCNE/C, 0, FA 

FON = (DEXP(-C*Z) + 0EXP(-0*Z) ) /( Z*Z + FA*FA) 

RETURN 

ENO 



OOUBLE PRECISION FUNCTION FCN(Z) 

FCN EVALUATES THE OIFFERENCE BETWEEN THE TRUE INTEGRANO ANO 
THE APPROXIMATE ONE FOR GIVEN Z 

IMPLICIT REAL*8 (A-H.O-Z) 

COMMON/ CON ST/ PHI, PH2.fr 1, FR2, ASCL , PI ,CA, Cl , AA 

COMMDN/GEZRO/=A,B,F,ICS,IOI 

COMMON/ RHO/N 

GO TO ( 1, 2, 1, 2 ) , ICS 

1 FR = FRl 
PH = PHI 
P = ASCL 
GO TO 3 

2 FR = FR2 
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PH = PH2 
P = O.lDl 

3 ARGF = FR*Z + AA=^=DLOG(Z) + PH 
ARGH = FA*Z + B 
S = 0. IDl - DEXP(-P=^^Z) 

ZN = 0.101/ (Z^=Z) 

NN = N + 2 
DO 6 I =1 ,NN 
6 ZN = ZN*Z 

ENV = S^S*S^DEXP(-Z)^ZN 
IF ( ICS .GT .2 ) GO TO 4 

FCN = ENV*(DSI N( ARGF) - FTDS IN ( ARGH) ) 
GO TO 5 

A FCN = ENV’^IDCOSC ARGF) - F*D S IN { ARGH ) ) 
5 RETURN 
END 
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SUBROUTINE D05 12 ( C , D , FCT , Y ) 



D0512 IS A 512-POINT G AUSS-L EGENDR E QUADRATURE WHICH EVAL- 
UATES THE BOUND-FREE INTEGRALS BETWEEN THE LIMITS C AND D. 
THE EXTERNALLY DEFINED FUNCTION FCT IS EITHER SIN OR COS. 
LOCATION OF THE SAMPLE POINTS AND THE CORRESPONDING WEIGHT 
FACTORS ARE READ INTO THE MAIN PROGRAM FROM AN EXTERNAL 
SOURCE . 

IMPLICIT REAL*=8 (A-H,0-Z) 

COMMON/ L51 2/ X( 256) , A(256) 

COMMON /CON ST/PHl ,PH2 , FR 1 , F R2 , A SC L , PI ,CA ,CT ,AA 

COMMON/GEZRO/F A, B, AMP , ICS , IQI 

COMMON /RHO/N 

GO, TO ( 1 ,2 , 1 ,2 ) , ICS 

1 FR = FRl 
PH = PHI 
P = ASCL 
GO TO 3 

2 FR = FR2 
PH = PH2 
P O.lDl 

3 Y = D.ODT 

CHANGE THE VARIABLE SO THAT THE INTEGRATION LIMITS BECCME 
(-1 . + 1 ) 

DPC = 0.5D0*(D + C) 

DMC = 0.5D0*(D - C) 

DO 4 I =1,256 

XDC = X ( I ) *dm: 

ZP = DPC + XDC 
ZM = DPC - XDC 
ZMN = 0. 1D1/(ZM*ZM) 

ZPN = 0. IDl/ (ZP* ZPO 
NN = N + 2 
DC 5 J=1,NN 
ZMN = ZMN«ZM 
5 ZPN = ZPN*ZP 

SP = O.lDl - DEXP(-P*ZP) 

SM = O.lDl - DEXP(-P*ZM) 
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ENVP = SP*SP-SP=!^DEXP(-ZP )*ZPN 
EtJVfl = SM^SM=5'5M^DEXP (-ZM )>-'ZMN 
ARGP = FR*ZP *- AA*DLCG(ZP) + PH 
ARGM = FR=fZM + AA=!=DLOG( ZI-1 ) + PH 
4 Y = Y + A( I )*( ENVP*FCT(ARGP) + ENVM--^FCT ( ARGM) ) 
Y = Y*DMC 
RETURN 
END 
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rHCsjrO 



SUBR3UT INE HnT(Sl,C 2 fSS,CC,I MAX, MAXI ) 



HINT EVALUATES THE BOUND-FREE INTEGRALS FOR THE CASE Z=1 
IMPLICIT REAL’f'8 (A-H,0-Z) 

EXTERNAL FACT 

COMMON/ CON ST/ PHI , PH2 , FRl , F R2 , AS CL , P I , C A , C I , A AA 

DIMENSION RHOl ( 41 , RH02 (4) , THTA Sl( 4) , THTACl (4) , 

1 THTAS2 (4 ) ,THTSC2 (4 ) 

DIMENSION Sin MAX) , C2( IMAX) , SS ( MA X I ) , CC ( MA XI) 

FI = FR1=^^=!'2 

F2 = FR2**2 

DO 1 K=i,4 

A = DF LOAT ( K - 1 ) 

RHOKK) = DS0RT((0.1D1 + A'*:<ASCL )^«2 + FI) 

RH02(K) = DSQRT((0.1D1 + A)*=f'2 + F2 ) 

THTASl(K) = DARS IN(FR1/RH01(K) ) 

THTACKK) = DARCOSI ( 0. IDl + A=^ A SC L ) / RH 01 ( K ) ) 

THTAS2(K) = OARS !N(FR2/RH02 (K) ) 

1 THTAC2(K) = DARCOS ( ( 0. ID 1 + A)/RHC2(K)) 

SKI) = DATAN(FRl) - 0 . 3 Dl=^ { DAT AN ( FRl / (0 . 1 D1 + ASCD) 

1 - DAT A\( FRl/ { O.lDl + 0 . 2C 1* A SC L ) ) ) 

2 - DATAN( FRl/ (D.lDl + J -a D1':ASCL ) ) 

C2(i) = 0. 15DK^OLOG( ( ( 0. ID 1 + ASCL)’f^*2 +• F1)/({0.1D1 

1 + 0. 2DK ASCL )^="2 + FD) + J . SDO^DLOG ( ( O . ID 1 

2 + 0. 3DK ASCL )^-'-2 + F1)/(0.1D1 + FD) 

SS(1) = DATAN(FR2) - 0 . 3 D1 ='' ( DAT AN ( FR2/ 0 * 2D 1 ) 

1 - DATAN( FR2/0.3D1 ) ) - DATAN ( FR2/0 . 4D1 ) 

CC(1) = 0. 15DDDL3G( (0.4D1 + F2)/(J.9D1 + F2)) + 

1 0. 5D0>"DL0G( ( 0. 16D2 + F2)/10.1D1 + F2 ) ) 

DO 2 J =2, I MAX 

I = J - 1 

FI = DFLOAT { I ) 

FCT = FACT(J - 2) 

SKJ) = FCT* ( DSI N( FDTHTASl ( 1 ) ) / RHOi (1 )^^>!=I 

- 0.3D1«{DSIN{FI*THTAS1(2) )/RH01(2)**I 

- DSIN( FDTHTASl (3 ) )/RH01( 3)^*1 ) 

- DSIN(Fr4THTASl(4) ) /RHOl (4)** I ) 

C2U) = FCT>^(DC0S(FI*THTAC1{1) )/RH01(l)**I 

- 0.3D1*{DCOS( FKTHTACK2) ) /RHOl (2)*# I 

- DCOS( FI*THTAC1( 3) )/RH01(3)**I ) 

- DCOS( F I*THTAC1 (4) ) / RHOl (4 )**I ) 

IF( J.GT.MAXI ) GO TO 2 



l6l 






SS(J) = 

- 0 . 



CC(J) = 
- 0 . 



2 CONTINUE 
RETURN 
END 



FCT=^= (DSIN( FI*THTAS2( 1 ) ) /RH02 (L 
3D1^(DSIN( FI=^THTAS2( 2) ) /RH02(2)*«I 

- DSIN(FI^THTAS2(3) ) / RH02 ( 3 ) 1 ) 

- DSIN(FI*THTAS2(4) J / RH02 ( 4 ) '-1= * I) 

FCT^’= ( DCOS ( F I=-=T HT AC2 ( 1 ) ) / RH02 ( 1 ) I 
3D1- (DCOS( FP"THTAC2(2) ) /RH02 (2)*=^ I 

- DCOS( FI-THTAC2I J) ) / RH02 ( 3 ) ) 

- DCOSIF I*THTAC2(4) ) / RH02 ( 4 ) I ) 
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If- <(■ -it- 



SUBROUTINE A 12 SUM ( N U, L MDA , AA ,VAL) 

A12SUM EVALUATES A 1 ( N U . L AMDA , 2 ) OR A2 ( NU , LAMDA , 2 ) USING II 
OR 13, RESPECTIVELY 

IMPLICIT REAL-i'S (A-H,0-Z) 

DI MENS ION VAL( 1 ) 

COMMON /TRM/ALFA , ET A 

COMMON/ATNO/Z 

ALFA2 = Z + ALFA 

FCTR = 0. 2D1* A L F A2 / ( 0. 2D 1=!'-Z + ALFA) 

NUl = NU + 1 

NL = NU + LMDA 

F = 0. IDl/DFLOAT ( NU) 

FCT = O.lDl/FCTR 
TVAL = O.ODO 
DO 10 1=1, NU 
F = F*DFL0AT(NU1 - I) 

FCT = FCT*FCTR 

10 TVAL = TVAL + FCT ^V AL ( NL-I ) *F 
AA = TVAL/ALFA2*"'( NL-2 ) 

RETURN 

END 



I- ^ 

y: yr * yt * 
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SU3R0UT INE PHASE ( AS , AC , BS , BC , H AX I , HS , HC , E ) 



PHASE COMPILES THE MATRIX ELEMENTS < S I ( H-E ) I CH I ( I) > AND 
<C 1 (H-E) i CHI ( I )> 

IMPLICIT REALMS (A-H,0-Z) 

COMMON /TRM/ALF A , ETA 

COMMON/ ATNO/Z 

COMMON/ I NDEX/\' I , N J , L I t L J , M I , MJ , I S 
DI ME NS I ON AS(1 ) , AC(1 I , BS (1 ) , BC ( 1 ) 

DIMENSION NU(26) ,LMDA( 26) , MU( 26) , COEF ( 26 ) 

FN = DFLOAT(NJ) 

FL = DFLOAT ( LJ ) 

FM = DFLCAT(MJ) 

SET THE VALUES OF NJ , LAMDA, MU AND THE COEFFICIENTS FOR 
EACH OF THE 26 TERMS OF THE MATRIX ELEMENT 



NU{ 1) = 


LJ + 1 


NU( 3) = 


NUI 1 ) 


NU( 4) = 


NUI 1) 


NU( 6) = 


NUIl) 


NUi e) = 


NUI 1) 


NU( 9) = 


NUI 1) 


NUI 19) = 


Null) 


LMDAI 12) 


= riui 1 ) 


LMDAI26) 


= NUIl) 


NU I 2 ) = 


LJ 


LMDAI 17) 


= NUI 2) 


NUI 21) = 


NU 12 ) 


NUI 5) = 


LJ + 2 


NUI20) = 


NUI5) 


NU I 2 2 ) = 


NUI 5) 


LMDAI 10) 


= NUI 5) 


LMDAIll ) 


= NUI 5) 


LMDAI 14) 


= NUI 5) 


LMDAI15) 


= NUI 5 ) 


LMDAI 16 ) 


= NUI 5) 


LMDAI 18) 


= NUI 5) 



LMDA(24) 
NU( 7) = 
NU { 1 0 ) = 
NU(12) = 
NU(13) = 

NU ( 1 5 ) = 
NU(17) = 
N'U (18) = 

NU(23) = 
LMDA( 3) 
LMDA(22 ) 
NU.(11) = 
NU(25) = 
LMDA( 8) 
NU(14) = 

NU(24) = 
NU(26) = 
LHDA( 1) 
LMDA( 2) 
LMDA( 5) 
LMDA( 6) 
LMOA( 7) 
LMDA( 9) 
LMDA( 20 ) 
NU(16) = 
LM0A( 4) 
LMDA( 19 ) 
LMDA( 21 ) 
LKDA(13 ) 
Lr'.DA(23) 
LMDA(25) 
MU( 1) = 
NU( 2) = 

MU( 3) = 
MU( 7) = 
HU ( 8 ) = 



= NU( 5 ) 

LJ - 1 
NJ + 1 
NU( 10 ) 

NU( 10 ) 

NU( 10 ) 

NU( 10 ) 

NU( 10 ) 

NU( 10 ) 

= NU( 10) 

= NU ( 1 0 ) 

NJ 

NUdl ) 

= NU( 11 ) 

NJ + 2 
NU ( 14 ) 

NU( 14) 

= NU( 14 ) 

= NU( 14) 

= NU( 14) 

= NU( 14 ) 

= NU( 14) 

= NU( 14 ) 

= r;U( 14 ) 

NJ - 1 
= NJ + 3 
= LMDA(4) 

= LM0A(4) 

= LJ + 3 
= LM0A( 13) 
= LMDA(13) 
MJ + 2 
MU( 1) 

MU ( 1 ) 

MU( 1) 

MU ( 1 ) 



MU( 10 ) 




HU( 1) • 


MU( 11) 


= 


MU( 1) 


MU (12) 




MU( 1 ) 


HU ( 1 6 ) 


= 


MU( 1) 


MU( 17) 


- 


MU ( 1 ) 


MU( 4) 


= 


MJ 


HU( 5) 


= 


HU(4) 


MU( 6) 


= 


MU ( 4 ) 


MU( 13) 




MU( 4) 


MU (14) 


= 


HU (4) 


MU( 15 ) 


- 


HU( 4) 


MU( 19) 


= 


MU( 4) 


MU (20) 


=r 


MU(4) 


MU(21) 


= 


MU( 4) 


MU (22) 


= 


MU ( 4 ) 


MU ( 23 ) 




MU( 4) 


MU( 24) 


= 


MU (4) 


HU (25 ) 


= 


MU (4 ) 


MU( 26) 


= 


MU( 4) 


MU( 9) 


= 


MJ + 1 


MU ( 18 ) 


= 


MU (9) 


COEF( 


1) 


= -0. ZEDO'" ALFA=^'v 


COEFdO ) 


= COEF( 1 ) 


COEF( 


2) 


= 0. 5D0^ALFA''M FL 


COEF ( 12 ) 


= COEF ( 2 ) 


COEF( 


3) 


= 0.5D0^ALFA* ( FN 


COEF( 1 1) 


= COEF ( 3) 


COEF( 


4 ) 


= 0.5 D0*ALFA*FM 


COEF( 


5) 


= C0EF(4) 


COEF ( 1 


3) 


= COEF ( 4 ) 


COEF( 14) 


= C0EF(4) 


COEF (19) 


= -COEF ( 4) 


COEF (20 ) 


= -CDEF( 4) 


C0EF(23) 


= -C0EF(4) 


COEF (24) 


= -COEF (4 ) 


COEF ( 


6 ) 


= -FM«(FM + FM + 



E 

. IDl ) - Z 
.101) - Z 



+ O.lDl) 
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COEF( 15) 
COEF( 7) 
COEF( 17) 
COEF( 8) 
COEF( 16 ) 
COEF( 9) 
C0EF(18) 
CQEF( 21) 
C0EF(26) 
COEF ( 22) 
COEF( 25) 
IF'( IS . EO 



= COEF ( 6 ) 

= -0. 5DO'-^^FL^^( FL 
= COEF( 7) 

= -0. FN-^ ( FN 
= COEF( 6) 

= 0.101 

= COEF( 9) 

= FM*FL 
= COEF (21 ) 

= FM^FN 
= C0EF(22) 

0) GO TO 10 



-{- O.lDl) 
+ 0 .101 ) 



CHAMGE THE SIGHS OF THE LAST 15 TERHS FOR THE TRIPLET CASE 
00 1 1=10,18 

1 COEF( I ) = -COEF( I ) 

DO 2 1=23,26 

2 COEF ( I ) = -COEF( I ) 

10 LLMAX = MAXI 

NHMAX = MAXI + 1 
HS = 0.0 00 
HC = 0.000 
00 20 K=l, 18 

IN = NU(K) + 2 + ^IN‘•'lAX^M LMOA( K ) + LLMAX*MU(K) - 1) 

HS = HS +- AS( I N)*COEF ( K) 

HC = HC + AC( I N)-COEF ( K ) 

20 CONTINUE 

DC 30 K=19,26 

IN = NU(K) + 2 + NNMAX*(LMDA(K) + LLHAX^KU(K) - 1) 

HS = HS + BS( I N)=!=COEF ( K) 

HC = HC + BC( IN)=^COEF(K ) 

30 CONTINUE 
RETURN 
END 
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location, but not the width, of the scattering resonances. Singlet and 
triplet s-wave phase shifts for e-H and e-He"^ scattering are compared with 
the results of other calculations and H* and He S state energy levels 
including the auto-ionizing levels are presented and compared with other 
calculations and with experiment. It is tentatively concluded that the 
Harris method does not work on systems whose long range potential is of 
the Coulomb type. 
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